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, with 

357 pixels in nJy (`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 identity = afwImage.PhotoCalib(1.0) 

454 if photoCalib is not None and photoCalib != identity: 

455 template.maskedImage = photoCalib.calibrateImage(template.maskedImage) 

456 photoCalib = identity 

457 template.setPhotoCalib(photoCalib) 

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

459 

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

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

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

463 template.getInfo().setCoaddInputs(coaddInputs) 

464 return pipeBase.Struct(template=template) 

465 

466 def _correctVariance(self, template, plateScaleFactor): 

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

468 

469 Compute a multiplicative correction factor for the template variance 

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

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

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

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

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

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

476 image pixel noise level would suggest. 

477 

478 Parameters 

479 ---------- 

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

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

482 plateScaleFactor : `float` or `None` 

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

484 

485 Raises 

486 ------ 

487 RuntimeError 

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

489 could not be reconstructed from the coadd inputs. 

490 """ 

491 scale = 1.0 

492 if self.config.doCorrectVariancePlateScale: 

493 if plateScaleFactor is None: 

494 raise RuntimeError( 

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

496 ) 

497 scale *= plateScaleFactor 

498 self.metadata["variancePlateScaleFactor"] = plateScaleFactor 

499 self.log.info( 

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

501 " coadd inputs.", 

502 plateScaleFactor, 

503 ) 

504 

505 if scale == 1.0: 

506 return 

507 template.variance.array *= scale 

508 self.metadata["templateVarianceCorrectionFactor"] = scale 

509 self.log.info( 

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

511 scale, 

512 ) 

513 

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

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

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

517 

518 Parameters 

519 ---------- 

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

521 `lsst.afw.image.Exposure` 

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

523 component is read; the pixels are left alone. 

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

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

526 bbox : `lsst.geom.Box2I` 

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

528 evaluate the science image pixel scale. 

529 

530 Returns 

531 ------- 

532 factor : `float` or `None` 

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

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

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

536 """ 

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

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

539 return None 

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

541 

542 for record in coaddInputs.ccds: 

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

544 # that gives a usable pixel area. 

545 originalWcs = record.getWcs() 

546 if originalWcs is None: 

547 continue 

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

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

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

551 factor = (pviScale/scienceScale)**2 

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

553 return factor 

554 return None 

555 

556 def checkHighVariance(self, template): 

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

558 

559 Parameters 

560 ---------- 

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

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

563 """ 

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

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

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

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

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

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

570 " pixels were unmasked for background estimation, but" 

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

572 else: 

573 varianceExposure = template.clone() 

574 varianceExposure.image.array = varianceExposure.variance.array 

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

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

577 highVariancePix = varianceBackground > threshold 

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

579 

580 @staticmethod 

581 def _checkInputs(dataIds, coaddExposures): 

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

583 the exposures all have the same photometric calibration. 

584 

585 Parameters 

586 ---------- 

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

588 Record of the tract and patch of each coaddExposure. 

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

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

591 `lsst.afw.image.Exposure` or 

592 `lsst.afw.image.Exposure`]] 

593 Coadds to be mosaicked. 

594 

595 Returns 

596 ------- 

597 band : `str` 

598 Filter band of all the input exposures. 

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

600 Photometric calibration of all of the input exposures. 

601 

602 Raises 

603 ------ 

604 RuntimeError 

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

606 all the same. 

607 """ 

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

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

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

611 band = bands.pop() 

612 photoCalibs = [ 

613 exposure.get(component="photoCalib") 

614 for exposures in coaddExposures.values() 

615 for exposure in exposures 

616 ] 

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

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

619 raise RuntimeError(msg) 

620 photoCalib = photoCalibs[0] 

621 return band, photoCalib 

622 

623 def _makeExposureCatalog(self, exposureRefs, dataIds): 

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

625 

626 Parameters 

627 ---------- 

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

629 `lsst.afw.image.Exposure`] 

630 Exposures to include in the catalog. 

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

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

633 "patch" entries. 

634 

635 Returns 

636 ------- 

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

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

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

640 Catalog of metadata for each exposure 

641 totalBox : `lsst.geom.Box2I` 

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

643 """ 

644 catalog = afwTable.ExposureCatalog(self.schema) 

645 catalog.reserve(len(exposureRefs)) 

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

647 images = {} 

648 totalBox = geom.Box2I() 

649 

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

651 images[dataId] = coadd.maskedImage 

652 bbox = coadd.getBBox() 

653 totalBox = totalBox.expandedTo(bbox) 

654 record = catalog.addNew() 

655 record.setPsf(coadd.psf) 

656 record.setWcs(coadd.wcs) 

657 record.setPhotoCalib(coadd.photoCalib) 

658 record.setBBox(bbox) 

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

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

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

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

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

664 record.set("weight", 1) 

665 

666 return images, catalog, totalBox 

667 

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

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

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

671 exposures. 

672 

673 Parameters 

674 ---------- 

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

676 `lsst.afw.image.Exposure`] 

677 Images to be merged into one larger bounding box. 

678 bbox : `lsst.geom.Box2I` 

679 Bounding box defining the image to merge into. 

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

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

682 

683 Returns 

684 ------- 

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

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

687 positions. 

688 count : `int` 

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

690 in the merged image. 

691 included : `list` [`int`] 

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

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

694 """ 

695 merged = afwImage.ExposureF(bbox, wcs) 

696 weights = afwImage.ImageF(bbox) 

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

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

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

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

701 clippedBox.clip(bbox) 

702 if clippedBox.area == 0: 

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

704 continue # nothing in this image overlaps the output 

705 maskedImage = maskedImage.subset(clippedBox) 

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

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

708 np.isfinite(maskedImage.variance.array) 

709 ) 

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

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

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

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

714 maskedImage.image.array[bad] = 0.0 

715 maskedImage.variance.array[bad] = 0.0 

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

717 maskedImage.mask.array[bad] = 0 

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

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

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

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

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

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

724 # Free memory before creating new large arrays 

725 del weight 

726 merged.maskedImage[clippedBox] += maskedImage 

727 included.append(i) 

728 

729 good = weights.array > 0 

730 

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

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

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

734 weights = weights.array[good] 

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

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

737 

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

739 

740 return merged, good.sum(), included 

741 

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

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

744 

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

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

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

748 went in to those pixels. 

749 

750 Parameters 

751 ---------- 

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

753 Generated template the PSF is for. 

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

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

756 of the input PSFs. 

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

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

759 

760 Returns 

761 ------- 

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

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

764 """ 

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

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

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

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

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

770 

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

772 coaddPsf = CoaddPsf( 

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

774 ) 

775 return coaddPsf 

776 

777 

778class GetDcrTemplateConnections( 

779 GetTemplateConnections, 

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

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

782): 

783 visitInfo = pipeBase.connectionTypes.Input( 

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

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

786 storageClass="VisitInfo", 

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

788 ) 

789 dcrCoadds = pipeBase.connectionTypes.Input( 

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

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

792 storageClass="ExposureF", 

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

794 multiple=True, 

795 deferLoad=True, 

796 ) 

797 

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

799 super().__init__(config=config) 

800 self.inputs.remove("coaddExposures") 

801 

802 

803class GetDcrTemplateConfig( 

804 GetTemplateConfig, pipelineConnections=GetDcrTemplateConnections 

805): 

806 numSubfilters = pexConfig.Field( 

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

808 dtype=int, 

809 default=3, 

810 ) 

811 effectiveWavelength = pexConfig.Field( 

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

813 optional=False, 

814 dtype=float, 

815 ) 

816 bandwidth = pexConfig.Field( 

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

818 optional=False, 

819 dtype=float, 

820 ) 

821 

822 def validate(self): 

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

824 raise ValueError( 

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

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

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

828 ) 

829 

830 

831class GetDcrTemplateTask(GetTemplateTask): 

832 ConfigClass = GetDcrTemplateConfig 

833 _DefaultName = "getDcrTemplate" 

834 

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

836 inputs = butlerQC.get(inputRefs) 

837 bbox = inputs.pop("bbox") 

838 wcs = inputs.pop("wcs") 

839 dcrCoaddExposureHandles = inputs.pop("dcrCoadds") 

840 skymap = inputs.pop("skyMap") 

841 visitInfo = inputs.pop("visitInfo") 

842 

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

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

845 

846 results = self.getExposures( 

847 dcrCoaddExposureHandles, bbox, skymap, wcs, visitInfo 

848 ) 

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

850 outputs = self.run( 

851 coaddExposureHandles=results.coaddExposures, 

852 bbox=bbox, 

853 wcs=wcs, 

854 dataIds=results.dataIds, 

855 physical_filter=physical_filter, 

856 ) 

857 butlerQC.put(outputs, outputRefs) 

858 

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

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

861 the detector. 

862 

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

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

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

866 

867 Parameters 

868 ---------- 

869 dcrCoaddExposureHandles : `list` \ 

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

871 `lsst.afw.image.Exposure`] 

872 Data references to exposures that might overlap the detector. 

873 bbox : `lsst.geom.Box2I` 

874 Template Bounding box of the detector geometry onto which to 

875 resample the coaddExposures. 

876 skymap : `lsst.skymap.SkyMap` 

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

878 template exposures. 

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

880 Template WCS onto which to resample the coaddExposures. 

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

882 Metadata for the science image. 

883 

884 Returns 

885 ------- 

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

887 A struct with attibutes: 

888 

889 ``coaddExposures`` 

890 Dict of coadd exposures that overlap the projected bbox, 

891 indexed on tract id 

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

893 ``dataIds`` 

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

895 projected bbox, indexed on tract id 

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

897 

898 Raises 

899 ------ 

900 pipeBase.NoWorkFound