Coverage for python/lsst/ip/diffim/getTemplate.py: 67%

1# This file is part of ip_diffim. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21import collections 

22import dataclasses 

23 

24import numpy as np 

25 

26import lsst.afw.image as afwImage 

27import lsst.geom as geom 

28import lsst.afw.geom as afwGeom 

29from lsst.afw.image import VisitInfo 

30import lsst.afw.table as afwTable 

31from lsst.afw.math._warper import computeWarpedBBox 

32import lsst.afw.math as afwMath 

33import lsst.pex.config as pexConfig 

34import lsst.pipe.base as pipeBase 

35 

36from lsst.skymap import BaseSkyMap 

37from lsst.ip.diffim.dcrModel import DcrModel 

38from lsst.meas.algorithms import CoaddPsf, CoaddPsfConfig, SubtractBackgroundTask, ScaleVarianceTask 

39from lsst.utils.timer import timeMethod 

40 

41__all__ = [ 

42 "GetTemplateTask", 

43 "GetTemplateConfig", 

44 "GetDcrTemplateTask", 

45 "GetDcrTemplateConfig", 

46] 

47 

48 

49class GetTemplateConnections( 

50 pipeBase.PipelineTaskConnections, 

51 dimensions=("instrument", "visit", "detector"), 

52 defaultTemplates={"coaddName": "goodSeeing", "warpTypeSuffix": "", "fakesType": ""}, 

53): 

54 bbox = pipeBase.connectionTypes.Input( 

55 doc="Bounding box of exposure to determine the geometry of the output template.", 

56 name="{fakesType}calexp.bbox", 

57 storageClass="Box2I", 

58 dimensions=("instrument", "visit", "detector"), 

59 ) 

60 wcs = pipeBase.connectionTypes.Input( 

61 doc="WCS of the exposure that we will construct the template for.", 

62 name="{fakesType}calexp.wcs", 

63 storageClass="Wcs", 

64 dimensions=("instrument", "visit", "detector"), 

65 ) 

66 skyMap = pipeBase.connectionTypes.Input( 

67 doc="Geometry of the tracts and patches that the coadds are defined on.", 

68 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

69 dimensions=("skymap",), 

70 storageClass="SkyMap", 

71 ) 

72 coaddExposures = pipeBase.connectionTypes.Input( 

73 doc="Coadds that may overlap the desired region, as possible inputs to the template." 

74 " Will be restricted to those that directly overlap the projected bounding box.", 

75 dimensions=("tract", "patch", "skymap", "band"), 

76 storageClass="ExposureF", 

77 name="{fakesType}{coaddName}Coadd{warpTypeSuffix}", 

78 multiple=True, 

79 deferLoad=True, 

80 deferGraphConstraint=True, 

81 ) 

82 

83 template = pipeBase.connectionTypes.Output( 

84 doc="Warped template, pixel matched to the bounding box and WCS.", 

85 dimensions=("instrument", "visit", "detector"), 

86 storageClass="ExposureF", 

87 name="{fakesType}{coaddName}Diff_templateExp{warpTypeSuffix}", 

88 ) 

89 

90 def __init__(self, *, config=None): 

91 super().__init__(config=config) 

92 if config.requireCoaddAtGraphBuild: 

93 self.coaddExposures = dataclasses.replace( 

94 self.coaddExposures, 

95 deferGraphConstraint=False, 

96 ) 

97 

98 

99class GetTemplateConfig( 

100 pipeBase.PipelineTaskConfig, pipelineConnections=GetTemplateConnections 

101): 

102 templateBorderSize = pexConfig.Field( 

103 dtype=int, 

104 default=20, 

105 doc="Number of pixels to grow the requested template image to account for warping", 

106 ) 

107 warp = pexConfig.ConfigField( 

108 dtype=afwMath.Warper.ConfigClass, 

109 doc="warper configuration", 

110 ) 

111 doCorrectVariancePlateScale = pexConfig.Field( 

112 dtype=bool, 

113 default=False, 

114 doc="Correct the template variance for the total change in pixel" 

115 " area between the images the coadds were built from and the science" 

116 " image? This matters when the two have different plate scales, such" 

117 " as DECam templates used for Rubin science images. The coadd pixel" 

118 " grid cancels out of the total, so the skymap does not enter. It is a" 

119 " single scalar factor, applied to the final template variance plane.", 

120 ) 

121 coaddPsf = pexConfig.ConfigField( 

122 doc="Configuration for CoaddPsf", 

123 dtype=CoaddPsfConfig, 

124 ) 

125 varianceBackground = pexConfig.ConfigurableField( 

126 target=SubtractBackgroundTask, 

127 doc="Task to estimate the background variance.", 

128 ) 

129 highVarianceThreshold = pexConfig.RangeField( 

130 dtype=float, 

131 default=4, 

132 min=1, 

133 doc="Set the HIGH_VARIANCE mask plane for regions with variance" 

134 " greater than the median by this factor.", 

135 ) 

136 highVarianceMaskFraction = pexConfig.Field( 

137 dtype=float, 

138 default=0.1, 

139 doc="Minimum fraction of unmasked pixels needed to set the" 

140 " HIGH_VARIANCE mask plane.", 

141 ) 

142 doScaleVariance = pexConfig.Field( 

143 dtype=bool, 

144 default=True, 

145 doc="Scale variance of the template image?" 

146 ) 

147 scaleVariance = pexConfig.ConfigurableField( 

148 target=ScaleVarianceTask, 

149 doc="Subtask to rescale the variance of the template to the statistically expected level." 

150 ) 

151 requireCoaddAtGraphBuild = pexConfig.Field( 

152 dtype=bool, 

153 default=False, 

154 doc="If True, include the coadd dataset existence in the" 

155 " initial butler query during QuantumGraph generation.", 

156 ) 

157 

158 def setDefaults(self): 

159 # Use a smaller cache: per SeparableKernel.computeCache, this should 

160 # give a warping error of a fraction of a count (these must match). 

161 self.warp.cacheSize = 100000 

162 self.coaddPsf.cacheSize = self.warp.cacheSize 

163 # The WCS for LSST should be smoothly varying, so we can use a longer 

164 # interpolation length for WCS evaluations. 

165 self.warp.interpLength = 100 

166 self.warp.warpingKernelName = "lanczos3" 

167 self.coaddPsf.warpingKernelName = self.warp.warpingKernelName 

168 

169 # Background subtraction of the variance plane 

170 self.varianceBackground.algorithm = "LINEAR" 

171 self.varianceBackground.binSize = 32 

172 self.varianceBackground.useApprox = False 

173 self.varianceBackground.statisticsProperty = "MEDIAN" 

174 self.varianceBackground.doFilterSuperPixels = True 

175 self.varianceBackground.ignoredPixelMask = ["BAD", 

176 "EDGE", 

177 "DETECTED", 

178 "DETECTED_NEGATIVE", 

179 "NO_DATA", 

180 ] 

181 

182 

183class GetTemplateTask(pipeBase.PipelineTask): 

184 ConfigClass = GetTemplateConfig 

185 _DefaultName = "getTemplate" 

186 

187 def __init__(self, *args, **kwargs): 

188 super().__init__(*args, **kwargs) 

189 if self.config.doScaleVariance: 

190 self.makeSubtask("scaleVariance") 

191 self.warper = afwMath.Warper.fromConfig(self.config.warp) 

192 self.schema = afwTable.ExposureTable.makeMinimalSchema() 

193 self.schema.addField( 

194 "tract", type=np.int32, doc="Which tract this exposure came from." 

195 ) 

196 self.schema.addField( 

197 "patch", 

198 type=np.int32, 

199 doc="Which patch in the tract this exposure came from.", 

200 ) 

201 self.schema.addField( 

202 "weight", 

203 type=float, 

204 doc="Weight for each exposure, used to make the CoaddPsf; should always be 1.", 

205 ) 

206 self.makeSubtask("varianceBackground") 

207 

208 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

209 inputs = butlerQC.get(inputRefs) 

210 bbox = inputs.pop("bbox") 

211 wcs = inputs.pop("wcs") 

212 coaddExposures = inputs.pop("coaddExposures") 

213 skymap = inputs.pop("skyMap") 

214 

215 # This should not happen with a properly configured execution context. 

216 assert not inputs, "runQuantum got more inputs than expected" 

217 

218 results = self.getExposures(coaddExposures, bbox, skymap, wcs) 

219 physical_filter = butlerQC.quantum.dataId["physical_filter"] 

220 outputs = self.run( 

221 coaddExposureHandles=results.coaddExposures, 

222 bbox=bbox, 

223 wcs=wcs, 

224 dataIds=results.dataIds, 

225 physical_filter=physical_filter, 

226 visit=outputRefs.template.dataId["visit"], 

227 ) 

228 butlerQC.put(outputs, outputRefs) 

229 

230 def getExposures(self, coaddExposureHandles, bbox, skymap, wcs): 

231 """Return a data structure containing the coadds that overlap the 

232 specified bbox projected onto the sky, and a corresponding data 

233 structure of their dataIds. 

234 These are the appropriate inputs to this task's `run` method. 

235 

236 The spatial index in the butler registry has generous padding and often 

237 supplies patches near, but not directly overlapping the desired region. 

238 This method filters the inputs so that `run` does not have to read in 

239 all possibly-matching coadd exposures. 

240 

241 Parameters 

242 ---------- 

243 coaddExposureHandles : `iterable` \ 

244 [`lsst.daf.butler.DeferredDatasetHandle` of \ 

245 `lsst.afw.image.Exposure`] 

246 Dataset handles to exposures that might overlap the desired 

247 region. 

248 bbox : `lsst.geom.Box2I` 

249 Template bounding box of the pixel geometry onto which the 

250 coaddExposures will be resampled. 

251 skymap : `lsst.skymap.SkyMap` 

252 Geometry of the tracts and patches the coadds are defined on. 

253 wcs : `lsst.afw.geom.SkyWcs` 

254 Template WCS onto which the coadds will be resampled. 

255 

256 Returns 

257 ------- 

258 result : `lsst.pipe.base.Struct` 

259 A struct with attributes: 

260 

261 ``coaddExposures`` 

262 Dict of coadd exposures that overlap the projected bbox, 

263 indexed on tract id 

264 (`dict` [`int`, `list` [`lsst.daf.butler.DeferredDatasetHandle` of 

265 `lsst.afw.image.Exposure`] ]). 

266 ``dataIds`` 

267 Dict of data IDs of the coadd exposures that overlap the 

268 projected bbox, indexed on tract id 

269 (`dict` [`int`, `list [`lsst.daf.butler.DataCoordinate`] ]). 

270 

271 Raises 

272 ------ 

273 NoWorkFound 

274 Raised if no patches overlap the input detector bbox, or the input 

275 WCS is None. 

276 """ 

277 if wcs is None: 

278 raise pipeBase.NoWorkFound( 

279 "WCS is None; cannot find overlapping exposures." 

280 ) 

281 

282 # Exposure's validPolygon would be more accurate 

283 detectorPolygon = geom.Box2D(bbox) 

284 detectorCorners = wcs.pixelToSky(detectorPolygon.getCorners()) 

285 overlappingArea = 0 

286 coaddExposures = collections.defaultdict(list) 

287 dataIds = collections.defaultdict(list) 

288 

289 for coaddRef in coaddExposureHandles: 

290 dataId = coaddRef.dataId 

291 patchWcs = skymap[dataId["tract"]].getWcs() 

292 patchBBox = skymap[dataId["tract"]][dataId["patch"]].getOuterBBox() 

293 patchPolygon = afwGeom.Polygon(geom.Box2D(patchBBox)) 

294 # Calculate detector/patch overlap in patch coordinates rather than 

295 # detector coordinates because the skymap's inverse mapping 

296 # (patchWcs.skyToPixel()) is more stable than the detector's for 

297 # arbitrary sky coordinates. 

298 detectorInPatchCoordinates = afwGeom.Polygon(patchWcs.skyToPixel(detectorCorners)) 

299 if patchPolygon.intersection(detectorInPatchCoordinates): 

300 overlappingArea += patchPolygon.intersectionSingle( 

301 detectorInPatchCoordinates 

302 ).calculateArea() 

303 self.log.info( 

304 "Using template input tract=%s, patch=%s", 

305 dataId["tract"], 

306 dataId["patch"], 

307 ) 

308 coaddExposures[dataId["tract"]].append(coaddRef) 

309 dataIds[dataId["tract"]].append(dataId) 

310 

311 if not overlappingArea: 

312 raise pipeBase.NoWorkFound("No patches overlap detector") 

313 

314 return pipeBase.Struct(coaddExposures=coaddExposures, dataIds=dataIds) 

315 

316 @timeMethod 

317 def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visit=None): 

318 """Warp coadds from multiple tracts and patches to form a template to 

319 subtract from a science image. 

320 

321 Tract and patch overlap regions are combined by a variance-weighted 

322 average, and the variance planes are combined with the same weights, 

323 not added in quadrature; the overlap regions are not statistically 

324 independent, because they're derived from the same original data. 

325 The PSF on the template is created by combining the CoaddPsf on each 

326 template image into a meta-CoaddPsf. 

327 

328 Parameters 

329 ---------- 

330 coaddExposureHandles : `dict` [`int`, `list` of \ 

331 [`lsst.daf.butler.DeferredDatasetHandle` of \ 

332 `lsst.afw.image.Exposure`]] 

333 Coadds to be mosaicked, indexed on tract id. 

334 bbox : `lsst.geom.Box2I` 

335 Template Bounding box of the detector geometry onto which to 

336 resample the ``coaddExposureHandles``. Modified in-place to include the 

337 template border. 

338 wcs : `lsst.afw.geom.SkyWcs` 

339 Template WCS onto which to resample the ``coaddExposureHandles``. 

340 dataIds : `dict` [`int`, `list` [`lsst.daf.butler.DataCoordinate`]] 

341 Record of the tract and patch of each coaddExposure, indexed on 

342 tract id. 

343 physical_filter : `str` 

344 Physical filter of the science image. 

345 visit : `int`, optional 

346 If supplied, over-write the visit ID in the template's visitInfo 

347 so that downstream source injection tasks can link the template and 

348 science image for the visit. 

349 

350 Returns 

351 ------- 

352 result : `lsst.pipe.base.Struct` 

353 A struct with attributes: 

354 

355 ``template`` 

356 A template coadd exposure assembled out of patches 

357 (`lsst.afw.image.ExposureF`). 

358 

359 Raises 

360 ------ 

361 NoWorkFound 

362 If no coadds are found with sufficient un-masked pixels. 

363 """ 

364 band, photoCalib = self._checkInputs(dataIds, coaddExposureHandles) 

365 

366 bbox.grow(self.config.templateBorderSize) 

367 

368 warped = {} 

369 catalogs = [] 

370 # Determine the ratio of the original pixel area to the pixel area of 

371 # the science image, if configured. This will only be different from 1 

372 # if the coadd was from a different instrument. The ratio can be 

373 # determined once from the components of one coadd exposure, so that no 

374 # pixels have to be read yet. Coadds comprising images from multiple 

375 # different instruments are not supported. 

376 plateScaleFactor = None 

377 if self.config.doCorrectVariancePlateScale and coaddExposureHandles: 

378 plateScaleFactor = self._plateScaleFactor( 

379 next(iter(coaddExposureHandles.values()))[0], wcs, bbox 

380 ) 

381 

382 for tract in coaddExposureHandles: 

383 maskedImages, catalog, totalBox = self._makeExposureCatalog( 

384 coaddExposureHandles[tract], dataIds[tract] 

385 ) 

386 warpedBox = computeWarpedBBox(catalog[0].wcs, bbox, wcs) 

387 warpedBox.grow(5) # to ensure we catch all relevant input pixels 

388 # Combine images from individual patches together. 

389 unwarped, count, included = self._merge( 

390 maskedImages, warpedBox, catalog[0].wcs 

391 ) 

392 # Delete `maskedImages` after combining into one large image to reduce peak memory use 

393 del maskedImages 

394 if count == 0: 

395 self.log.info( 

396 "No valid pixels from coadd patches in tract %s; not including in output.", 

397 tract, 

398 ) 

399 continue 

400 warpedBox.clip(totalBox) 

401 potentialInput = self.warper.warpExposure( 

402 wcs, unwarped.subset(warpedBox), destBBox=bbox 

403 ) 

404 # Delete the single large `unwarped` image after warping to reduce peak memory use 

405 del unwarped 

406 if np.all( 

407 potentialInput.mask.array 

408 & potentialInput.mask.getPlaneBitMask("NO_DATA") 

409 ): 

410 self.log.info( 

411 "No overlap from coadd patches in tract %s; not including in output.", 

412 tract, 

413 ) 

414 continue 

415 

416 # Trim the exposure catalog to just the patches that were used. 

417 tempCatalog = afwTable.ExposureCatalog(self.schema) 

418 tempCatalog.reserve(len(included)) 

419 for i in included: 

420 tempCatalog.append(catalog[i]) 

421 catalogs.append(tempCatalog) 

422 warped[tract] = potentialInput.maskedImage 

423 

424 if len(warped) == 0: 

425 raise pipeBase.NoWorkFound("No patches found to overlap science exposure.") 

426 

427 # At this point, all entries will be valid, so we can ignore included. 

428 template, count, _ = self._merge(warped, bbox, wcs) 

429 if count == 0: 429 ↛ 430line 429 didn't jump to line 430 because the condition on line 429 was never true

430 raise pipeBase.NoWorkFound("No valid pixels in warped template.") 

431 

432 if self.config.doScaleVariance: 

433 # Scale the variance of the template image before subtraction, if 

434 # needed. Note that the science variance is scaled 

435 # independently in ``AlardLuptonSubtractTask``. 

436 varianceFactor = self.scaleVariance.run(template.maskedImage) 

437 self.log.info("Template variance scaling factor: %.2f", varianceFactor) 

438 self.metadata["scaleTemplateVarianceFactor"] = varianceFactor 

439 

440 self._correctVariance(template, plateScaleFactor) 

441 

442 # Make a single catalog containing all the inputs that were accepted. 

443 catalog = afwTable.ExposureCatalog(self.schema) 

444 catalog.reserve(sum([len(c) for c in catalogs])) 

445 for c in catalogs: 

446 catalog.extend(c) 

447 

448 # Set a mask plane for any regions with exceptionally high variance. 

449 self.checkHighVariance(template) 

450 if visit is not None: 450 ↛ 451line 450 didn't jump to line 451 because the condition on line 450 was never true

451 template.getInfo().setVisitInfo(VisitInfo(id=visit)) 

452 template.setFilter(afwImage.FilterLabel(band, physical_filter)) 

453 template.setPhotoCalib(photoCalib) 

454 template.setPsf(self._makePsf(template, catalog, wcs)) 

455 

456 # Record the input coadd patches as the template's coadd inputs. 

457 coaddInputs = afwImage.CoaddInputs(afwTable.ExposureTable.makeMinimalSchema(), self.schema) 

458 coaddInputs.ccds.extend(catalog, deep=True) 

459 template.getInfo().setCoaddInputs(coaddInputs) 

460 return pipeBase.Struct(template=template) 

461 

462 def _correctVariance(self, template, plateScaleFactor): 

463 """Correct the template variance plane for the effects of warping. 

464 

465 Compute a multiplicative correction factor for the template variance 

466 if the plate scale of the coadd's constituent images is different than 

467 the science image the template is being constructed for. This should 

468 only be necessary if the coadd images were from a different instrument 

469 than the science image. If the plate scale of the science instrument is 

470 smaller than the plate scale of the coadd instrument, then the template 

471 pixels will be correlated and the true variance will be higher than the 

472 image pixel noise level would suggest. 

473 

474 Parameters 

475 ---------- 

476 template : `lsst.afw.image.Exposure` 

477 Assembled template; its variance plane is modified in place. 

478 plateScaleFactor : `float` or `None` 

479 Correction for the change in pixel area, from `_plateScaleFactor`. 

480 

481 Raises 

482 ------ 

483 RuntimeError 

484 Raised if ``doCorrectVariancePlateScale`` is set but the factor 

485 could not be reconstructed from the coadd inputs. 

486 """ 

487 scale = 1.0 

488 if self.config.doCorrectVariancePlateScale: 

489 if plateScaleFactor is None: 

490 raise RuntimeError( 

491 "doCorrectVariancePlateScale is set but no usable coaddInputs were found." 

492 ) 

493 scale *= plateScaleFactor 

494 self.metadata["variancePlateScaleFactor"] = plateScaleFactor 

495 self.log.info( 

496 "Applying a plate scale variance factor of %.4f, reconstructed from the" 

497 " coadd inputs.", 

498 plateScaleFactor, 

499 ) 

500 

501 if scale == 1.0: 

502 return 

503 template.variance.array *= scale 

504 self.metadata["templateVarianceCorrectionFactor"] = scale 

505 self.log.info( 

506 "Corrected the template variance plane for warping by a factor of %.4f.", 

507 scale, 

508 ) 

509 

510 def _plateScaleFactor(self, coaddHandle, wcs, bbox): 

511 """Compute the variance correction for the total change in pixel area 

512 between the coadd's constituent images and the science image. 

513 

514 Parameters 

515 ---------- 

516 coaddHandle : `lsst.daf.butler.DeferredDatasetHandle` of \ 

517 `lsst.afw.image.Exposure` 

518 Handle to one of the input coadd patches. Only its ``coaddInputs`` 

519 component is read; the pixels are left alone. 

520 wcs : `lsst.afw.geom.SkyWcs` 

521 WCS of the science image the template is being built for. 

522 bbox : `lsst.geom.Box2I` 

523 Bounding box of the template, used only to choose where to 

524 evaluate the science image pixel scale. 

525 

526 Returns 

527 ------- 

528 factor : `float` or `None` 

529 The correction factor equivalent to the ratio of the area of a pixel 

530 from the coadd's instrument to the science instrument, or `None` if 

531 this patch carries no input record from which it could be computed. 

532 """ 

533 coaddInputs = coaddHandle.get(component="coaddInputs") 

534 if coaddInputs is None: 534 ↛ 535line 534 didn't jump to line 535 because the condition on line 534 was never true

535 return None 

536 scienceScale = wcs.getPixelScale(geom.Box2D(bbox).getCenter()).asArcseconds() 

537 

538 for record in coaddInputs.ccds: 

539 # Iterate through the records, and use the first one with a WCS 

540 # that gives a usable pixel area. 

541 originalWcs = record.getWcs() 

542 if originalWcs is None: 

543 continue 

544 center = geom.Box2D(record.getBBox()).getCenter() 

545 pviScale = originalWcs.getPixelScale(center).asArcseconds() 

546 # The area of one original pixel, expressed in science pixels. 

547 factor = (pviScale/scienceScale)**2 

548 if np.isfinite(factor) and factor > 0: 548 ↛ 538line 548 didn't jump to line 538 because the condition on line 548 was always true

549 return factor 

550 return None 

551 

552 def checkHighVariance(self, template): 

553 """Set a mask plane for regions with unusually high variance. 

554 

555 Parameters 

556 ---------- 

557 template : `lsst.afw.image.Exposure` 

558 The warped template exposure, which will be modified in place. 

559 """ 

560 highVarianceMaskPlaneBit = template.mask.addMaskPlane("HIGH_VARIANCE") 

561 ignoredPixelBits = template.mask.getPlaneBitMask(self.varianceBackground.config.ignoredPixelMask) 

562 goodMask = (template.mask.array & ignoredPixelBits) == 0 

563 goodFraction = np.count_nonzero(goodMask)/template.mask.array.size 

564 if goodFraction < self.config.highVarianceMaskFraction: 564 ↛ 565line 564 didn't jump to line 565 because the condition on line 564 was never true

565 self.log.info("Not setting HIGH_VARIANCE mask plane, only %2.1f%% of" 

566 " pixels were unmasked for background estimation, but" 

567 " %2.1f%% are required", 100*goodFraction, 100*self.config.highVarianceMaskFraction) 

568 else: 

569 varianceExposure = template.clone() 

570 varianceExposure.image.array = varianceExposure.variance.array 

571 varianceBackground = self.varianceBackground.run(varianceExposure).background.getImage().array 

572 threshold = self.config.highVarianceThreshold*np.nanmedian(varianceBackground) 

573 highVariancePix = varianceBackground > threshold 

574 template.mask.array[highVariancePix] |= 2**highVarianceMaskPlaneBit 

575 

576 @staticmethod 

577 def _checkInputs(dataIds, coaddExposures): 

578 """Check that the all the dataIds are from the same band and that 

579 the exposures all have the same photometric calibration. 

580 

581 Parameters 

582 ---------- 

583 dataIds : `dict` [`int`, `list` [`lsst.daf.butler.DataCoordinate`]] 

584 Record of the tract and patch of each coaddExposure. 

585 coaddExposures : `dict` [`int`, `list` of \ 

586 [`lsst.daf.butler.DeferredDatasetHandle` of \ 

587 `lsst.afw.image.Exposure` or 

588 `lsst.afw.image.Exposure`]] 

589 Coadds to be mosaicked. 

590 

591 Returns 

592 ------- 

593 band : `str` 

594 Filter band of all the input exposures. 

595 photoCalib : `lsst.afw.image.PhotoCalib` 

596 Photometric calibration of all of the input exposures. 

597 

598 Raises 

599 ------ 

600 RuntimeError 

601 Raised if the bands or calibrations of the input exposures are not 

602 all the same. 

603 """ 

604 bands = set(dataId["band"] for tract in dataIds for dataId in dataIds[tract]) 

605 if len(bands) > 1: 605 ↛ 606line 605 didn't jump to line 606 because the condition on line 605 was never true

606 raise RuntimeError(f"GetTemplateTask called with multiple bands: {bands}") 

607 band = bands.pop() 

608 photoCalibs = [ 

609 exposure.get(component="photoCalib") 

610 for exposures in coaddExposures.values() 

611 for exposure in exposures 

612 ] 

613 if not all([photoCalibs[0] == x for x in photoCalibs]): 613 ↛ 614line 613 didn't jump to line 614 because the condition on line 613 was never true

614 msg = f"GetTemplateTask called with exposures with different photoCalibs: {photoCalibs}" 

615 raise RuntimeError(msg) 

616 photoCalib = photoCalibs[0] 

617 return band, photoCalib 

618 

619 def _makeExposureCatalog(self, exposureRefs, dataIds): 

620 """Make an exposure catalog for one tract. 

621 

622 Parameters 

623 ---------- 

624 exposureRefs : `list` of [`lsst.daf.butler.DeferredDatasetHandle` of \ 

625 `lsst.afw.image.Exposure`] 

626 Exposures to include in the catalog. 

627 dataIds : `list` [`lsst.daf.butler.DataCoordinate`] 

628 Data ids of each of the included exposures; must have "tract" and 

629 "patch" entries. 

630 

631 Returns 

632 ------- 

633 images : `dict` [`lsst.afw.image.MaskedImage`] 

634 MaskedImages of each of the input exposures, for warping. 

635 catalog : `lsst.afw.table.ExposureCatalog` 

636 Catalog of metadata for each exposure 

637 totalBox : `lsst.geom.Box2I` 

638 The union of the bounding boxes of all the input exposures. 

639 """ 

640 catalog = afwTable.ExposureCatalog(self.schema) 

641 catalog.reserve(len(exposureRefs)) 

642 exposures = (exposureRef.get() for exposureRef in exposureRefs) 

643 images = {} 

644 totalBox = geom.Box2I() 

645 

646 for coadd, dataId in zip(exposures, dataIds): 

647 images[dataId] = coadd.maskedImage 

648 bbox = coadd.getBBox() 

649 totalBox = totalBox.expandedTo(bbox) 

650 record = catalog.addNew() 

651 record.setPsf(coadd.psf) 

652 record.setWcs(coadd.wcs) 

653 record.setPhotoCalib(coadd.photoCalib) 

654 record.setBBox(bbox) 

655 record.setValidPolygon(afwGeom.Polygon(geom.Box2D(bbox).getCorners())) 

656 record.set("tract", dataId["tract"]) 

657 record.set("patch", dataId["patch"]) 

658 # Weight is used by CoaddPsf, but the PSFs from overlapping patches 

659 # should be very similar, so this value mostly shouldn't matter. 

660 record.set("weight", 1) 

661 

662 return images, catalog, totalBox 

663 

664 def _merge(self, maskedImages, bbox, wcs): 

665 """Merge the images that came from one tract into one larger image, 

666 ignoring NaN pixels and non-finite variance pixels from individual 

667 exposures. 

668 

669 Parameters 

670 ---------- 

671 maskedImages : `dict` [`lsst.afw.image.MaskedImage` or 

672 `lsst.afw.image.Exposure`] 

673 Images to be merged into one larger bounding box. 

674 bbox : `lsst.geom.Box2I` 

675 Bounding box defining the image to merge into. 

676 wcs : `lsst.afw.geom.SkyWcs` 

677 WCS of all of the input images to set on the output image. 

678 

679 Returns 

680 ------- 

681 merged : `lsst.afw.image.MaskedImage` 

682 Merged image with all of the inputs at their respective bbox 

683 positions. 

684 count : `int` 

685 Count of the number of good pixels (those with positive weights) 

686 in the merged image. 

687 included : `list` [`int`] 

688 List of indexes of patches that were included in the merged 

689 result, to be used to trim the exposure catalog. 

690 """ 

691 merged = afwImage.ExposureF(bbox, wcs) 

692 weights = afwImage.ImageF(bbox) 

693 included = [] # which patches were included in the result 

694 for i, (dataId, maskedImage) in enumerate(maskedImages.items()): 

695 # Only merge into the trimmed box, to save memory 

696 clippedBox = geom.Box2I(maskedImage.getBBox()) 

697 clippedBox.clip(bbox) 

698 if clippedBox.area == 0: 

699 self.log.debug("%s does not overlap template region.", dataId) 

700 continue # nothing in this image overlaps the output 

701 maskedImage = maskedImage.subset(clippedBox) 

702 # Catch both zero-value and NaN variance plane pixels 

703 good = (maskedImage.variance.array > 0) & ( 

704 np.isfinite(maskedImage.variance.array) 

705 ) 

706 weight = maskedImage.variance.array[good] ** (-0.5) 

707 bad = np.isnan(maskedImage.image.array) | ~good 

708 # Note that modifying the patch MaskedImage in place is fine; 

709 # we're throwing it away at the end anyway. 

710 maskedImage.image.array[bad] = 0.0 

711 maskedImage.variance.array[bad] = 0.0 

712 # Reset mask, too, since these pixels don't contribute to sum. 

713 maskedImage.mask.array[bad] = 0 

714 # Cannot use `merged.maskedImage *= weight` because that operator 

715 # multiplies the variance by the weight twice; in this case 

716 # `weight` are the exact values we want to scale by. 

717 maskedImage.image.array[good] *= weight 

718 maskedImage.variance.array[good] *= weight 

719 weights[clippedBox].array[good] += weight 

720 # Free memory before creating new large arrays 

721 del weight 

722 merged.maskedImage[clippedBox] += maskedImage 

723 included.append(i) 

724 

725 good = weights.array > 0 

726 

727 # Cannot use `merged.maskedImage /= weights` because that 

728 # operator divides the variance by the weight twice; in this case 

729 # `weights` are the exact values we want to scale by. 

730 weights = weights.array[good] 

731 merged.image.array[good] /= weights 

732 merged.variance.array[good] /= weights 

733 

734 merged.mask.array[~good] |= merged.mask.getPlaneBitMask("NO_DATA") 

735 

736 return merged, good.sum(), included 

737 

738 def _makePsf(self, template, catalog, wcs): 

739 """Return a PSF containing the PSF at each of the input regions. 

740 

741 Note that although this includes all the exposures from the catalog, 

742 the PSF knows which part of the template the inputs came from, so when 

743 evaluated at a given position it will not include inputs that never 

744 went in to those pixels. 

745 

746 Parameters 

747 ---------- 

748 template : `lsst.afw.image.Exposure` 

749 Generated template the PSF is for. 

750 catalog : `lsst.afw.table.ExposureCatalog` 

751 Catalog of exposures that went into the template that contains all 

752 of the input PSFs. 

753 wcs : `lsst.afw.geom.SkyWcs` 

754 WCS of the template, to warp the PSFs to. 

755 

756 Returns 

757 ------- 

758 coaddPsf : `lsst.meas.algorithms.CoaddPsf` 

759 The meta-psf constructed from all of the input catalogs. 

760 """ 

761 # CoaddPsf centroid not only must overlap image, but must overlap the 

762 # part of image with data. Use centroid of region with data. 

763 boolmask = template.mask.array & template.mask.getPlaneBitMask("NO_DATA") == 0 

764 maskx = afwImage.makeMaskFromArray(boolmask.astype(afwImage.MaskPixel)) 

765 centerCoord = afwGeom.SpanSet.fromMask(maskx, 1).computeCentroid() 

766 

767 ctrl = self.config.coaddPsf.makeControl() 

768 coaddPsf = CoaddPsf( 

769 catalog, wcs, centerCoord, ctrl.warpingKernelName, ctrl.cacheSize 

770 ) 

771 return coaddPsf 

772 

773 

774class GetDcrTemplateConnections( 

775 GetTemplateConnections, 

776 dimensions=("instrument", "visit", "detector"), 

777 defaultTemplates={"coaddName": "dcr", "warpTypeSuffix": "", "fakesType": ""}, 

778): 

779 visitInfo = pipeBase.connectionTypes.Input( 

780 doc="VisitInfo of calexp used to determine observing conditions.", 

781 name="{fakesType}calexp.visitInfo", 

782 storageClass="VisitInfo", 

783 dimensions=("instrument", "visit", "detector"), 

784 ) 

785 dcrCoadds = pipeBase.connectionTypes.Input( 

786 doc="Input DCR template to match and subtract from the exposure", 

787 name="{fakesType}dcrCoadd{warpTypeSuffix}", 

788 storageClass="ExposureF", 

789 dimensions=("tract", "patch", "skymap", "band", "subfilter"), 

790 multiple=True, 

791 deferLoad=True, 

792 ) 

793 

794 def __init__(self, *, config=None): 

795 super().__init__(config=config) 

796 self.inputs.remove("coaddExposures") 

797 

798 

799class GetDcrTemplateConfig( 

800 GetTemplateConfig, pipelineConnections=GetDcrTemplateConnections 

801): 

802 numSubfilters = pexConfig.Field( 

803 doc="Number of subfilters in the DcrCoadd.", 

804 dtype=int, 

805 default=3, 

806 ) 

807 effectiveWavelength = pexConfig.Field( 

808 doc="Effective wavelength of the filter in nm.", 

809 optional=False, 

810 dtype=float, 

811 ) 

812 bandwidth = pexConfig.Field( 

813 doc="Bandwidth of the physical filter.", 

814 optional=False, 

815 dtype=float, 

816 ) 

817 

818 def validate(self): 

819 if self.effectiveWavelength is None or self.bandwidth is None: 

820 raise ValueError( 

821 "The effective wavelength and bandwidth of the physical filter " 

822 "must be set in the getTemplate config for DCR coadds. " 

823 "Required until transmission curves are used in DM-13668." 

824 ) 

825 

826 

827class GetDcrTemplateTask(GetTemplateTask): 

828 ConfigClass = GetDcrTemplateConfig 

829 _DefaultName = "getDcrTemplate" 

830 

831 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

832 inputs = butlerQC.get(inputRefs) 

833 bbox = inputs.pop("bbox") 

834 wcs = inputs.pop("wcs") 

835 dcrCoaddExposureHandles = inputs.pop("dcrCoadds") 

836 skymap = inputs.pop("skyMap") 

837 visitInfo = inputs.pop("visitInfo") 

838 

839 # This should not happen with a properly configured execution context. 

840 assert not inputs, "runQuantum got more inputs than expected" 

841 

842 results = self.getExposures( 

843 dcrCoaddExposureHandles, bbox, skymap, wcs, visitInfo 

844 ) 

845 physical_filter = butlerQC.quantum.dataId["physical_filter"] 

846 outputs = self.run( 

847 coaddExposureHandles=results.coaddExposures, 

848 bbox=bbox, 

849 wcs=wcs, 

850 dataIds=results.dataIds, 

851 physical_filter=physical_filter, 

852 ) 

853 butlerQC.put(outputs, outputRefs) 

854 

855 def getExposures(self, dcrCoaddExposureHandles, bbox, skymap, wcs, visitInfo): 

856 """Return lists of coadds and their corresponding dataIds that overlap 

857 the detector. 

858 

859 The spatial index in the registry has generous padding and often 

860 supplies patches near, but not directly overlapping the detector. 

861 Filters inputs so that we don't have to read in all input coadds. 

862 

863 Parameters 

864 ---------- 

865 dcrCoaddExposureHandles : `list` \ 

866 [`lsst.daf.butler.DeferredDatasetHandle` of \ 

867 `lsst.afw.image.Exposure`] 

868 Data references to exposures that might overlap the detector. 

869 bbox : `lsst.geom.Box2I` 

870 Template Bounding box of the detector geometry onto which to 

871 resample the coaddExposures. 

872 skymap : `lsst.skymap.SkyMap` 

873 Input definition of geometry/bbox and projection/wcs for 

874 template exposures. 

875 wcs : `lsst.afw.geom.SkyWcs` 

876 Template WCS onto which to resample the coaddExposures. 

877 visitInfo : `lsst.afw.image.VisitInfo` 

878 Metadata for the science image. 

879 

880 Returns 

881 ------- 

882 result : `lsst.pipe.base.Struct` 

883 A struct with attibutes: 

884 

885 ``coaddExposures`` 

886 Dict of coadd exposures that overlap the projected bbox, 

887 indexed on tract id 

888 (`dict` [`int`, `list` [`lsst.afw.image.Exposure`] ]). 

889 ``dataIds`` 

890 Dict of data IDs of the coadd exposures that overlap the 

891 projected bbox, indexed on tract id 

892 (`dict` [`int`, `list [`lsst.daf.butler.DataCoordinate`] ]). 

893 

894 Raises 

895 ------ 

896 pipeBase.NoWorkFound 

897 Raised if no patches overlatp the input detector bbox. 

898 """ 

899 # Check that the patches actually overlap the detector 

900 # Exposure's validPolygon would be more accurate 

901 if wcs is None: 

902 raise pipeBase.NoWorkFound("Exposure has no WCS; cannot create a template.") 

903