Loading [MathJax]/extensions/tex2jax.js
LSST Applications 28.0.0,g1653933729+a8ce1bb630,g1a997c3884+a8ce1bb630,g28da252d5a+5bd70b7e6d,g2bbee38e9b+638fca75ac,g2bc492864f+638fca75ac,g3156d2b45e+07302053f8,g347aa1857d+638fca75ac,g35bb328faa+a8ce1bb630,g3a166c0a6a+638fca75ac,g3e281a1b8c+7bbb0b2507,g4005a62e65+17cd334064,g414038480c+5b5cd4fff3,g41af890bb2+4ffae9de63,g4e1a3235cc+0f1912dca3,g6249c6f860+3c3976f90c,g80478fca09+46aba80bd6,g82479be7b0+77990446f6,g858d7b2824+78ba4d1ce1,g89c8672015+f667a5183b,g9125e01d80+a8ce1bb630,ga5288a1d22+2a6264e9ca,gae0086650b+a8ce1bb630,gb58c049af0+d64f4d3760,gc22bb204ba+78ba4d1ce1,gc28159a63d+638fca75ac,gcf0d15dbbd+32ddb6096f,gd6b7c0dfd1+3e339405e9,gda3e153d99+78ba4d1ce1,gda6a2b7d83+32ddb6096f,gdaeeff99f8+1711a396fd,gdd5a9049c5+b18c39e5e3,ge2409df99d+a5e4577cdc,ge33fd446bb+78ba4d1ce1,ge79ae78c31+638fca75ac,gf0baf85859+64e8883e75,gf5289d68f6+e1b046a8d7,gfa443fc69c+91d9ed1ecf,gfda6b12a05+8419469a56
LSST Data Management Base Package
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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
24
25import lsst.afw.image as afwImage
26import lsst.geom as geom
27import lsst.afw.geom as afwGeom
28import lsst.afw.table as afwTable
29import lsst.afw.math as afwMath
30import lsst.pex.config as pexConfig
31import lsst.pipe.base as pipeBase
32from lsst.skymap import BaseSkyMap
33from lsst.ip.diffim.dcrModel import DcrModel
34from lsst.meas.algorithms import CoaddPsf, CoaddPsfConfig
35from lsst.utils.timer import timeMethod
36
37__all__ = ["GetTemplateTask", "GetTemplateConfig",
38 "GetDcrTemplateTask", "GetDcrTemplateConfig"]
39
40
41class GetTemplateConnections(pipeBase.PipelineTaskConnections,
42 dimensions=("instrument", "visit", "detector", "skymap"),
43 defaultTemplates={"coaddName": "goodSeeing",
44 "warpTypeSuffix": "",
45 "fakesType": ""}):
46 bbox = pipeBase.connectionTypes.Input(
47 doc="BBoxes of calexp used determine geometry of output template",
48 name="{fakesType}calexp.bbox",
49 storageClass="Box2I",
50 dimensions=("instrument", "visit", "detector"),
51 )
52 wcs = pipeBase.connectionTypes.Input(
53 doc="WCS of the calexp that we want to fetch the template for",
54 name="{fakesType}calexp.wcs",
55 storageClass="Wcs",
56 dimensions=("instrument", "visit", "detector"),
57 )
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",
63 )
64 # TODO DM-31292: Add option to use global external wcs from jointcal
65 # Needed for DRP HSC
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}",
71 multiple=True,
72 deferLoad=True
73 )
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}",
79 )
80
81
82class GetTemplateConfig(pipeBase.PipelineTaskConfig,
83 pipelineConnections=GetTemplateConnections):
84 templateBorderSize = pexConfig.Field(
85 dtype=int,
86 default=20,
87 doc="Number of pixels to grow the requested template image to account for warping"
88 )
89 warp = pexConfig.ConfigField(
90 dtype=afwMath.Warper.ConfigClass,
91 doc="warper configuration",
92 )
93 coaddPsf = pexConfig.ConfigField(
94 doc="Configuration for CoaddPsf",
95 dtype=CoaddPsfConfig,
96 )
97
98 def setDefaults(self):
99 self.warp.warpingKernelName = 'lanczos5'
100 self.coaddPsf.warpingKernelName = 'lanczos5'
101
102
103class GetTemplateTask(pipeBase.PipelineTask):
104 ConfigClass = GetTemplateConfig
105 _DefaultName = "getTemplate"
106
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.')
115
116 def runQuantum(self, butlerQC, inputRefs, outputRefs):
117 inputs = butlerQC.get(inputRefs)
118 results = self.getOverlappingExposures(inputs)
119 del inputs["skyMap"] # Only needed for the above.
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)
125
126 def getOverlappingExposures(self, inputs):
127 """Return lists of coadds and their corresponding dataIds that overlap
128 the detector.
129
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.
133
134 Parameters
135 ----------
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
146 template exposures.
147 - wcs : `lsst.afw.geom.SkyWcs`
148 Template WCS onto which to resample the coaddExposures.
149
150 Returns
151 -------
152 result : `lsst.pipe.base.Struct`
153 A struct with attributes:
154
155 ``coaddExposures``
156 Dict of Coadd exposures that overlap the detector, indexed on
157 tract id (`dict` [`int`, `list` [`lsst.afw.image.Exposure`] ]).
158 ``dataIds``
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`] ]).
162
163 Raises
164 ------
165 NoWorkFound
166 Raised if no patches overlap the input detector bbox.
167 """
168 # Check that the patches actually overlap the detector
169 # Exposure's validPolygon would be more accurate
170 detectorPolygon = geom.Box2D(inputs['bbox'])
171 overlappingArea = 0
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:
181 patchPolygon = afwGeom.Polygon(inputsWcs.skyToPixel(patchCorners))
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)
188 else:
189 self.log.warning("Exposure %s has no WCS, so cannot include it in the template.",
190 coaddRef)
191
192 if not overlappingArea:
193 raise pipeBase.NoWorkFound('No patches overlap detector')
194
195 return pipeBase.Struct(coaddExposures=coaddExposures,
196 dataIds=dataIds)
197
198 @timeMethod
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.
202
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.
209
210 Parameters
211 ----------
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
217 template border.
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
222 tract id.
223 physical_filter : `str`
224 Physical filter of the science image.
225
226 Returns
227 -------
228 result : `lsst.pipe.base.Struct`
229 A struct with attributes:
230
231 ``template``
232 A template coadd exposure assembled out of patches
233 (`lsst.afw.image.ExposureF`).
234
235 Raises
236 ------
237 NoWorkFound
238 If no coadds are found with sufficient un-masked pixels.
239 """
240 band, photoCalib = self._checkInputs(dataIds, coaddExposures)
241
242 bbox.grow(self.config.templateBorderSize)
243
244 warped = {}
245 catalogs = []
246 for tract in coaddExposures:
247 maskedImages, catalog, totalBox = self._makeExposureCatalog(coaddExposures[tract],
248 dataIds[tract])
249 # Combine images from individual patches together.
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)
254 continue
255 catalogs.append(catalog)
256 warped[tract] = potentialInput
257 warped[tract].setWcs(wcs)
258
259 if len(warped) == 0:
260 raise pipeBase.NoWorkFound("No patches found to overlap science exposure.")
261 template = self._merge([x.maskedImage for x in warped.values()], bbox, wcs)
262
263 # Make a single catalog containing all the inputs that were accepted.
264 catalog = afwTable.ExposureCatalog(self.schema)
265 catalog.reserve(sum([len(c) for c in catalogs]))
266 for c in catalogs:
267 catalog.extend(c)
268
269 template.setPsf(self._makePsf(template, catalog, wcs))
270 template.setFilter(afwImage.FilterLabel(band, physical_filter))
271 template.setPhotoCalib(photoCalib)
272 return pipeBase.Struct(template=template)
273
274 @staticmethod
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.
278
279 Parameters
280 ----------
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.
285
286 Returns
287 -------
288 band : `str`
289 Filter band of all the input exposures.
290 photoCalib : `lsst.afw.image.PhotoCalib`
291 Photometric calibration of all of the input exposures.
292
293 Raises
294 ------
295 RuntimeError
296 Raised if the bands or calibrations of the input exposures are not
297 all the same.
298 """
299 bands = set(dataId["band"] for tract in dataIds for dataId in dataIds[tract])
300 if len(bands) > 1:
301 raise RuntimeError(f"GetTemplateTask called with multiple bands: {bands}")
302 band = bands.pop()
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
309
310 def _makeExposureCatalog(self, exposures, dataIds):
311 """Make an exposure catalog for one tract.
312
313 Parameters
314 ----------
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
319 "patch" entries.
320
321 Returns
322 -------
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.
329 """
330 catalog = afwTable.ExposureCatalog(self.schema)
331 catalog.reserve(len(exposures))
332 images = [exposure.maskedImage for exposure in exposures]
333 totalBox = geom.Box2I()
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())
341 record.setValidPolygon(afwGeom.Polygon(geom.Box2D(coadd.getBBox()).getCorners()))
342 record.set("tract", dataId["tract"])
343 record.set("patch", dataId["patch"])
344 # Weight is used by CoaddPsf, but the PSFs from overlapping patches
345 # should be very similar, so this value mostly shouldn't matter.
346 record.set("weight", 1)
347
348 return images, catalog, totalBox
349
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
353 exposures.
354
355 Parameters
356 ----------
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.
363
364 Returns
365 -------
366 merged : `lsst.afw.image.MaskedImage`
367 Merged image with all of the inputs at their respective bbox
368 positions.
369 """
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)
375 # Note that modifying the patch MaskedImage in place is fine;
376 # we're throwing it away at the end anyway.
377 maskedImage.image.array[bad] = 0.0
378 maskedImage.variance.array[bad] = 0.0
379 # Reset mask, too, since these pixels don't contribute to sum.
380 maskedImage.mask.array[bad] = 0
381 # Cannot use `merged.maskedImage *= weight` because that operator
382 # multiplies the variance by the weight twice; in this case
383 # `weight` are the exact values we want to scale by.
384 maskedImage.image *= weight
385 maskedImage.variance *= weight
386 merged.maskedImage[maskedImage.getBBox()] += maskedImage
387 # Clear the NaNs to ensure that areas missing from this input are
388 # masked with NO_DATA after the loop.
389 weight.array[np.isnan(weight.array)] = 0
390 weights[maskedImage.getBBox()] += weight
391 # Cannot use `merged.maskedImage /= weights` because that operator
392 # divides the variance by the weight twice; in this case `weights` are
393 # the exact values we want to scale by.
394 merged.image /= weights
395 merged.variance /= weights
396 merged.mask.array |= merged.mask.getPlaneBitMask("NO_DATA") * (weights.array == 0)
397
398 return merged
399
400 def _makePsf(self, template, catalog, wcs):
401 """Return a PSF containing the PSF at each of the input regions.
402
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.
407
408 Parameters
409 ----------
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
414 of the input PSFs.
415 wcs : `lsst.afw.geom.SkyWcs`
416 WCS of the template, to warp the PSFs to.
417
418 Returns
419 -------
420 coaddPsf : `lsst.meas.algorithms.CoaddPsf`
421 The meta-psf constructed from all of the input catalogs.
422 """
423 # CoaddPsf centroid not only must overlap image, but must overlap the
424 # part of image with data. Use centroid of region with data.
425 boolmask = template.mask.array & template.mask.getPlaneBitMask('NO_DATA') == 0
426 maskx = afwImage.makeMaskFromArray(boolmask.astype(afwImage.MaskPixel))
427 centerCoord = afwGeom.SpanSet.fromMask(maskx, 1).computeCentroid()
428
429 ctrl = self.config.coaddPsf.makeControl()
430 coaddPsf = CoaddPsf(catalog, wcs, centerCoord, ctrl.warpingKernelName, ctrl.cacheSize)
431 return coaddPsf
432
433
435 dimensions=("instrument", "visit", "detector", "skymap"),
436 defaultTemplates={"coaddName": "dcr",
437 "warpTypeSuffix": "",
438 "fakesType": ""}):
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"),
444 )
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"),
450 multiple=True,
451 deferLoad=True
452 )
453
454 def __init__(self, *, config=None):
455 super().__init__(config=config)
456 self.inputs.remove("coaddExposures")
457
458
459class GetDcrTemplateConfig(GetTemplateConfig,
460 pipelineConnections=GetDcrTemplateConnections):
461 numSubfilters = pexConfig.Field(
462 doc="Number of subfilters in the DcrCoadd.",
463 dtype=int,
464 default=3,
465 )
466 effectiveWavelength = pexConfig.Field(
467 doc="Effective wavelength of the filter.",
468 optional=False,
469 dtype=float,
470 )
471 bandwidth = pexConfig.Field(
472 doc="Bandwidth of the physical filter.",
473 optional=False,
474 dtype=float,
475 )
476
477 def validate(self):
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.")
482
483
484class GetDcrTemplateTask(GetTemplateTask):
485 ConfigClass = GetDcrTemplateConfig
486 _DefaultName = "getDcrTemplate"
487
488 def getOverlappingExposures(self, inputs):
489 """Return lists of coadds and their corresponding dataIds that overlap
490 the detector.
491
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.
495
496 Parameters
497 ----------
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
508 template exposures.
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.
513
514 Returns
515 -------
516 result : `lsst.pipe.base.Struct`
517 A struct with attibutes:
518
519 ``coaddExposures``
520 Coadd exposures that overlap the detector (`list`
521 [`lsst.afw.image.Exposure`]).
522 ``dataIds``
523 Data IDs of the coadd exposures that overlap the detector
524 (`list` [`lsst.daf.butler.DataCoordinate`]).
525
526 Raises
527 ------
528 NoWorkFound
529 Raised if no patches overlatp the input detector bbox.
530 """
531 # Check that the patches actually overlap the detector
532 # Exposure's validPolygon would be more accurate
533 detectorPolygon = geom.Box2D(inputs["bbox"])
534 overlappingArea = 0
535 coaddExposureRefList = []
536 dataIds = []
537 patchList = dict()
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())
543 patchPolygon = afwGeom.Polygon(inputs["wcs"].skyToPixel(patchCorners))
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'])
551 else:
552 patchList[dataId['tract']] = [dataId['patch'], ]
553 dataIds.append(dataId)
554
555 if not overlappingArea:
556 raise pipeBase.NoWorkFound('No patches overlap detector')
557
558 self.checkPatchList(patchList)
559
560 coaddExposures = self.getDcrModel(patchList, inputs['dcrCoadds'], inputs['visitInfo'])
561 return pipeBase.Struct(coaddExposures=coaddExposures,
562 dataIds=dataIds)
563
564 def checkPatchList(self, patchList):
565 """Check that all of the DcrModel subfilters are present for each
566 patch.
567
568 Parameters
569 ----------
570 patchList : `dict`
571 Dict of the patches containing valid data for each tract.
572
573 Raises
574 ------
575 RuntimeError
576 If the number of exposures found for a patch does not match the
577 number of subfilters.
578 """
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)
584
585 def getDcrModel(self, patchList, coaddRefs, visitInfo):
586 """Build DCR-matched coadds from a list of exposure references.
587
588 Parameters
589 ----------
590 patchList : `dict`
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.
597
598 Returns
599 -------
600 coaddExposures : `list` [`lsst.afw.image.Exposure`]
601 Coadd exposures that overlap the detector.
602 """
603 coaddExposures = []
604 for tract in patchList:
605 for patch in set(patchList[tract]):
606 coaddRefList = [coaddRef for coaddRef in coaddRefs
607 if _selectDataRef(coaddRef, tract, patch)]
608
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
615
616
617def _selectDataRef(coaddRef, tract, patch):
618 condition = (coaddRef.dataId['tract'] == tract) & (coaddRef.dataId['patch'] == patch)
619 return condition
Cartesian polygons.
Definition polygon.py:29
A group of labels for a filter in an exposure or coadd.
Definition FilterLabel.h:58
Custom catalog class for ExposureRecord/Table.
Definition fwd.h:67
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)
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)