Coverage for python / lsst / drp / tasks / dcr_assemble_coadd.py: 13%

1# This file is part of drp_tasks. 

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/>. 

21 

22__all__ = ["DcrAssembleCoaddConnections", "DcrAssembleCoaddTask", "DcrAssembleCoaddConfig"] 

23 

24from math import ceil 

25 

26import numpy as np 

27from scipy import ndimage 

28 

29import lsst.afw.image as afwImage 

30import lsst.afw.table as afwTable 

31import lsst.coadd.utils as coaddUtils 

32import lsst.geom as geom 

33import lsst.meas.algorithms as measAlg 

34import lsst.pex.config as pexConfig 

35import lsst.pipe.base as pipeBase 

36import lsst.utils as utils 

37from lsst.ip.diffim.dcrModel import DcrModel, applyDcr, calculateDcr 

38from lsst.meas.base import SingleFrameMeasurementTask 

39from lsst.pipe.tasks.coaddBase import makeSkyInfo, subBBoxIter 

40from lsst.pipe.tasks.measurePsf import MeasurePsfTask 

41from lsst.utils.timer import timeMethod 

42 

43from .assemble_coadd import ( 

44 CompareWarpAssembleCoaddConfig, 

45 CompareWarpAssembleCoaddConnections, 

46 CompareWarpAssembleCoaddTask, 

47) 

48 

49 

50class DcrAssembleCoaddConnections( 

51 CompareWarpAssembleCoaddConnections, 

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

53 defaultTemplates={ 

54 "inputWarpName": "deep", 

55 "inputCoaddName": "deep", 

56 "outputCoaddName": "dcr", 

57 "warpType": "direct", 

58 "warpTypeSuffix": "", 

59 "fakesType": "", 

60 }, 

61): 

62 inputWarps = pipeBase.connectionTypes.Input( 

63 doc=( 

64 "Input list of warps to be assembled i.e. stacked." 

65 "Note that this will often be different than the inputCoaddName." 

66 "WarpType (e.g. direct, psfMatched) is controlled by the warpType config parameter" 

67 ), 

68 name="{inputWarpName}Coadd_{warpType}Warp", 

69 storageClass="ExposureF", 

70 dimensions=("tract", "patch", "skymap", "visit", "instrument"), 

71 deferLoad=True, 

72 multiple=True, 

73 ) 

74 templateExposure = pipeBase.connectionTypes.Input( 

75 doc="Input coadded exposure, produced by previous call to AssembleCoadd", 

76 name="{fakesType}{inputCoaddName}Coadd{warpTypeSuffix}", 

77 storageClass="ExposureF", 

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

79 ) 

80 dcrCoadds = pipeBase.connectionTypes.Output( 

81 doc="Output coadded exposure, produced by stacking input warps", 

82 name="{fakesType}{outputCoaddName}Coadd{warpTypeSuffix}", 

83 storageClass="ExposureF", 

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

85 multiple=True, 

86 ) 

87 dcrNImages = pipeBase.connectionTypes.Output( 

88 doc="Output image of number of input images per pixel", 

89 name="{outputCoaddName}Coadd_nImage", 

90 storageClass="ImageU", 

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

92 multiple=True, 

93 ) 

94 

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

96 super().__init__(config=config) 

97 if not config.doWrite: 

98 self.outputs.remove("dcrCoadds") 

99 if not config.doNImage: 

100 self.outputs.remove("dcrNImages") 

101 # Remove outputs inherited from ``AssembleCoaddConnections`` that are 

102 # not used. 

103 self.outputs.remove("coaddExposure") 

104 self.outputs.remove("nImage") 

105 self.inputs.remove("psfMatchedWarps") 

106 

107 

108class DcrAssembleCoaddConfig(CompareWarpAssembleCoaddConfig, pipelineConnections=DcrAssembleCoaddConnections): 

109 dcrNumSubfilters = pexConfig.Field( 

110 dtype=int, 

111 doc="Number of sub-filters to forward model chromatic effects to fit the supplied exposures.", 

112 default=3, 

113 ) 

114 maxNumIter = pexConfig.Field( 

115 dtype=int, 

116 optional=True, 

117 doc="Maximum number of iterations of forward modeling.", 

118 default=None, 

119 ) 

120 minNumIter = pexConfig.Field( 

121 dtype=int, 

122 optional=True, 

123 doc="Minimum number of iterations of forward modeling.", 

124 default=None, 

125 ) 

126 convergenceThreshold = pexConfig.Field( 

127 dtype=float, 

128 doc="Target relative change in convergence between iterations of forward modeling.", 

129 default=0.001, 

130 ) 

131 useConvergence = pexConfig.Field( 

132 dtype=bool, 

133 doc="Use convergence test as a forward modeling end condition?" 

134 "If not set, skips calculating convergence and runs for ``maxNumIter`` iterations", 

135 default=True, 

136 ) 

137 baseGain = pexConfig.Field( 

138 dtype=float, 

139 optional=True, 

140 doc="Relative weight to give the new solution vs. the last solution when updating the model." 

141 "A value of 1.0 gives equal weight to both solutions." 

142 "Small values imply slower convergence of the solution, but can " 

143 "help prevent overshooting and failures in the fit." 

144 "If ``baseGain`` is None, a conservative gain " 

145 "will be calculated from the number of subfilters. ", 

146 default=None, 

147 ) 

148 useProgressiveGain = pexConfig.Field( 

149 dtype=bool, 

150 doc="Use a gain that slowly increases above ``baseGain`` to accelerate convergence? " 

151 "When calculating the next gain, we use up to 5 previous gains and convergence values." 

152 "Can be set to False to force the model to change at the rate of ``baseGain``. ", 

153 default=True, 

154 ) 

155 doAirmassWeight = pexConfig.Field( 

156 dtype=bool, 

157 doc="Weight exposures by airmass? Useful if there are relatively few high-airmass observations.", 

158 default=False, 

159 ) 

160 modelWeightsWidth = pexConfig.Field( 

161 dtype=float, 

162 doc="Width of the region around detected sources to include in the DcrModel.", 

163 default=3, 

164 ) 

165 useModelWeights = pexConfig.Field( 

166 dtype=bool, 

167 doc="Width of the region around detected sources to include in the DcrModel.", 

168 default=True, 

169 ) 

170 splitSubfilters = pexConfig.Field( 

171 dtype=bool, 

172 doc="Calculate DCR for two evenly-spaced wavelengths in each subfilter. Instead of at the midpoint", 

173 default=True, 

174 ) 

175 splitThreshold = pexConfig.Field( 

176 dtype=float, 

177 doc="Minimum DCR difference within a subfilter to use ``splitSubfilters``, in pixels." 

178 "Set to 0 to always split the subfilters.", 

179 default=0.1, 

180 ) 

181 regularizeModelIterations = pexConfig.Field( 

182 dtype=float, 

183 doc="Maximum relative change of the model allowed between iterations. Set to zero to disable.", 

184 default=2.0, 

185 ) 

186 regularizeModelFrequency = pexConfig.Field( 

187 dtype=float, 

188 doc="Maximum relative change of the model allowed between subfilters. Set to zero to disable.", 

189 default=4.0, 

190 ) 

191 convergenceMaskPlanes = pexConfig.ListField( 

192 dtype=str, default=["DETECTED"], doc="Mask planes to use to calculate convergence." 

193 ) 

194 regularizationWidth = pexConfig.Field( 

195 dtype=int, default=2, doc="Minimum radius of a region to include in regularization, in pixels." 

196 ) 

197 imageInterpOrder = pexConfig.Field( 

198 dtype=int, 

199 doc="The order of the spline interpolation used to shift the image plane.", 

200 default=1, 

201 ) 

202 accelerateModel = pexConfig.Field( 

203 dtype=float, 

204 doc="Factor to amplify the differences between model planes by to speed convergence.", 

205 default=3, 

206 ) 

207 doCalculatePsf = pexConfig.Field( 

208 dtype=bool, 

209 doc="Set to detect stars and recalculate the PSF from the final coadd." 

210 "Otherwise the PSF is estimated from a selection of the best input exposures", 

211 default=False, 

212 ) 

213 detectPsfSources = pexConfig.ConfigurableField( 

214 target=measAlg.SourceDetectionTask, 

215 doc="Task to detect sources for PSF measurement, if ``doCalculatePsf`` is set.", 

216 ) 

217 measurePsfSources = pexConfig.ConfigurableField( 

218 target=SingleFrameMeasurementTask, 

219 doc="Task to measure sources for PSF measurement, if ``doCalculatePsf`` is set.", 

220 ) 

221 measurePsf = pexConfig.ConfigurableField( 

222 target=MeasurePsfTask, 

223 doc="Task to measure the PSF of the coadd, if ``doCalculatePsf`` is set.", 

224 ) 

225 effectiveWavelength = pexConfig.Field( 

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

227 "Required if transmission curves aren't used." 

228 "Support for using transmission curves is to be added in DM-13668.", 

229 dtype=float, 

230 ) 

231 bandwidth = pexConfig.Field( 

232 doc="Bandwidth of the physical filter, in nm." 

233 "Required if transmission curves aren't used." 

234 "Support for using transmission curves is to be added in DM-13668.", 

235 dtype=float, 

236 ) 

237 

238 def setDefaults(self): 

239 CompareWarpAssembleCoaddConfig.setDefaults(self) 

240 self.doWriteArtifactMasks = False 

241 self.doNImage = True 

242 self.assembleStaticSkyModel.warpType = self.warpType 

243 # The coadd and nImage files will be overwritten by this 

244 # Task, so don't write them the first time. 

245 self.assembleStaticSkyModel.doNImage = False 

246 self.assembleStaticSkyModel.doWrite = False 

247 self.detectPsfSources.returnOriginalFootprints = False 

248 self.detectPsfSources.thresholdPolarity = "positive" 

249 # Only use bright sources for PSF measurement 

250 self.detectPsfSources.thresholdValue = 50 

251 self.detectPsfSources.nSigmaToGrow = 2 

252 # A valid star for PSF measurement should at least fill 5x5 pixels 

253 self.detectPsfSources.minPixels = 25 

254 # Use the variance plane to calculate signal to noise 

255 self.detectPsfSources.thresholdType = "pixel_stdev" 

256 # Ensure psf candidate size is as large as piff psf size. 

257 if ( 

258 self.doCalculatePsf 

259 and self.measurePsf.psfDeterminer.name == "piff" 

260 and self.psfDeterminer["piff"].kernelSize > self.makePsfCandidates.kernelSize 

261 ): 

262 self.makePsfCandidates.kernelSize = self.psfDeterminer["piff"].kernelSize 

263 

264 

265class DcrAssembleCoaddTask(CompareWarpAssembleCoaddTask): 

266 """Assemble DCR coadded images from a set of warps. 

267 

268 Attributes 

269 ---------- 

270 bufferSize : `int` 

271 The number of pixels to grow each subregion by to allow for DCR. 

272 

273 Notes 

274 ----- 

275 As with AssembleCoaddTask, we want to assemble a coadded image from a set 

276 of Warps (also called coadded temporary exposures), including the effects 

277 of Differential Chromatic Refraction (DCR). 

278 For full details of the mathematics and algorithm, please see 

279 DMTN-037: DCR-matched template generation (https://dmtn-037.lsst.io). 

280 

281 This Task produces a DCR-corrected Coadd, as well as a dcrCoadd 

282 for each subfilter used in the iterative calculation. 

283 It begins by dividing the bandpass-defining filter into N equal bandwidth 

284 "subfilters", and divides the flux in each pixel from an initial coadd 

285 equally into each as a "dcrModel". Because the airmass and parallactic 

286 angle of each individual exposure is known, we can calculate the shift 

287 relative to the center of the band in each subfilter due to DCR. For each 

288 exposure we apply this shift as a linear transformation to the dcrModels 

289 and stack the results to produce a DCR-matched exposure. The matched 

290 exposures are subtracted from the input exposures to produce a set of 

291 residual images, and these residuals are reverse shifted for each 

292 exposures' subfilters and stacked. The shifted and stacked residuals are 

293 added to the dcrModels to produce a new estimate of the flux in each pixel 

294 within each subfilter. The dcrModels are solved for iteratively, which 

295 continues until the solution from a new iteration improves by less than 

296 a set percentage, or a maximum number of iterations is reached. 

297 Two forms of regularization are employed to reduce unphysical results. 

298 First, the new solution is averaged with the solution from the previous 

299 iteration, which mitigates oscillating solutions where the model 

300 overshoots with alternating very high and low values. 

301 Second, a common degeneracy when the data have a limited range of airmass 

302 or parallactic angle values is for one subfilter to be fit with very low or 

303 negative values, while another subfilter is fit with very high values. This 

304 typically appears in the form of holes next to sources in one subfilter, 

305 and corresponding extended wings in another. Because each subfilter has 

306 a narrow bandwidth we assume that physical sources that are above the noise 

307 level will not vary in flux by more than a factor of `frequencyClampFactor` 

308 between subfilters, and pixels that have flux deviations larger than that 

309 factor will have the excess flux distributed evenly among all subfilters. 

310 If `splitSubfilters` is set, then each subfilter will be further sub- 

311 divided during the forward modeling step (only). This approximates using 

312 a higher number of subfilters that may be necessary for high airmass 

313 observations, but does not increase the number of free parameters in the 

314 fit. This is needed when there are high airmass observations which would 

315 otherwise have significant DCR even within a subfilter. Because calculating 

316 the shifted images takes most of the time, splitting the subfilters is 

317 turned off by way of the `splitThreshold` option for low-airmass 

318 observations that do not suffer from DCR within a subfilter. 

319 """ 

320 

321 ConfigClass = DcrAssembleCoaddConfig 

322 _DefaultName = "dcrAssembleCoadd" 

323 

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

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

326 if self.config.doCalculatePsf: 

327 self.schema = afwTable.SourceTable.makeMinimalSchema() 

328 self.makeSubtask("detectPsfSources", schema=self.schema) 

329 self.makeSubtask("measurePsfSources", schema=self.schema) 

330 self.makeSubtask("measurePsf", schema=self.schema) 

331 

332 @utils.inheritDoc(pipeBase.PipelineTask) 

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

334 # Docstring to be formatted with info from PipelineTask.runQuantum 

335 """ 

336 Notes 

337 ----- 

338 Assemble a coadd from a set of Warps. 

339 """ 

340 inputData = butlerQC.get(inputRefs) 

341 

342 # Construct skyInfo expected by run 

343 # Do not remove skyMap from inputData in case _makeSupplementaryData 

344 # needs it. 

345 skyMap = inputData["skyMap"] 

346 outputDataId = butlerQC.quantum.dataId 

347 

348 inputData["skyInfo"] = makeSkyInfo( 

349 skyMap, tractId=outputDataId["tract"], patchId=outputDataId["patch"] 

350 ) 

351 

352 # Construct list of input Deferred Datasets 

353 warpRefList = inputData["inputWarps"] 

354 

355 inputs = self.prepareInputs(warpRefList, inputData["skyInfo"].bbox) 

356 self.log.info("Found %d %s", len(inputs.warpRefList), self.getTempExpDatasetName(self.warpType)) 

357 if len(inputs.warpRefList) == 0: 

358 self.log.warning("No coadd temporary exposures found") 

359 return 

360 

361 supplementaryData = self._makeSupplementaryData(butlerQC, inputRefs, outputRefs) 

362 retStruct = self.run( 

363 inputData["skyInfo"], 

364 warpRefList=inputs.warpRefList, 

365 imageScalerList=inputs.imageScalerList, 

366 weightList=inputs.weightList, 

367 supplementaryData=supplementaryData, 

368 ) 

369 

370 inputData.setdefault("brightObjectMask", None) 

371 for subfilter in range(self.config.dcrNumSubfilters): 

372 # Use the PSF of the stacked dcrModel, and do not recalculate the 

373 # PSF for each subfilter 

374 retStruct.dcrCoadds[subfilter].setPsf(retStruct.coaddExposure.getPsf()) 

375 self.processResults(retStruct.dcrCoadds[subfilter], inputData["brightObjectMask"], outputDataId) 

376 

377 if self.config.doWrite: 

378 butlerQC.put(retStruct, outputRefs) 

379 return retStruct 

380 

381 def _makeSupplementaryData(self, butlerQC, inputRefs, outputRefs): 

382 """Load the previously-generated template coadd. 

383 

384 Returns 

385 ------- 

386 templateCoadd : `lsst.pipe.base.Struct` 

387 Results as a struct with attributes: 

388 

389 ``templateCoadd`` 

390 Coadded exposure (`lsst.afw.image.ExposureF`). 

391 """ 

392 templateCoadd = butlerQC.get(inputRefs.templateExposure) 

393 

394 return pipeBase.Struct(templateCoadd=templateCoadd) 

395 

396 def measureCoaddPsf(self, coaddExposure): 

397 """Detect sources on the coadd exposure and measure the final PSF. 

398 

399 Parameters 

400 ---------- 

401 coaddExposure : `lsst.afw.image.Exposure` 

402 The final coadded exposure. 

403 """ 

404 table = afwTable.SourceTable.make(self.schema) 

405 detResults = self.detectPsfSources.run(table, coaddExposure, clearMask=False) 

406 coaddSources = detResults.sources 

407 self.measurePsfSources.run(measCat=coaddSources, exposure=coaddExposure) 

408 # Measure the PSF on the stacked subfilter coadds if possible. 

409 # We should already have a decent estimate of the coadd PSF, however, 

410 # so in case of any errors simply log them as a warning and use the 

411 # default PSF. 

412 try: 

413 psfResults = self.measurePsf.run(coaddExposure, coaddSources) 

414 except Exception as e: 

415 self.log.warning("Unable to calculate PSF, using default coadd PSF: %s", e) 

416 else: 

417 coaddExposure.setPsf(psfResults.psf) 

418 

419 def prepareDcrInputs(self, templateCoadd, warpRefList, weightList): 

420 """Prepare the DCR coadd by iterating through the visitInfo of the 

421 input warps. 

422 

423 Sets the property ``bufferSize``. 

424 

425 Parameters 

426 ---------- 

427 templateCoadd : `lsst.afw.image.ExposureF` 

428 The initial coadd exposure before accounting for DCR. 

429 warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` 

430 The data references to the input warped exposures. 

431 weightList : `list` of `float` 

432 The weight to give each input exposure in the coadd. 

433 Will be modified in place if ``doAirmassWeight`` is set. 

434 

435 Returns 

436 ------- 

437 dcrModels : `lsst.pipe.tasks.DcrModel` 

438 Best fit model of the true sky after correcting chromatic effects. 

439 

440 Raises 

441 ------ 

442 NotImplementedError 

443 If ``lambdaMin`` is missing from the Mapper class of the obs 

444 package being used. 

445 """ 

446 sigma2fwhm = 2.0 * np.sqrt(2.0 * np.log(2.0)) 

447 filterLabel = templateCoadd.getFilter() 

448 dcrShifts = [] 

449 airmassDict = {} 

450 angleDict = {} 

451 psfSizeDict = {} 

452 for visitNum, warpExpRef in enumerate(warpRefList): 

453 visitInfo = warpExpRef.get(component="visitInfo") 

454 psf = warpExpRef.get(component="psf") 

455 visit = warpExpRef.dataId["visit"] 

456 # Just need a rough estimate; average positions are fine 

457 psfAvgPos = psf.getAveragePosition() 

458 psfSize = psf.computeShape(psfAvgPos).getDeterminantRadius() * sigma2fwhm 

459 airmass = visitInfo.getBoresightAirmass() 

460 parallacticAngle = visitInfo.getBoresightParAngle().asDegrees() 

461 airmassDict[str(visit)] = airmass 

462 angleDict[str(visit)] = parallacticAngle 

463 psfSizeDict[str(visit)] = psfSize 

464 if self.config.doAirmassWeight: 

465 weightList[visitNum] *= airmass 

466 dcrShifts.append( 

467 np.max( 

468 np.abs( 

469 calculateDcr( 

470 visitInfo, 

471 templateCoadd.getWcs(), 

472 self.config.effectiveWavelength, 

473 self.config.bandwidth, 

474 self.config.dcrNumSubfilters, 

475 ) 

476 ) 

477 ) 

478 ) 

479 self.log.info("Selected airmasses:\n%s", airmassDict) 

480 self.log.info("Selected parallactic angles:\n%s", angleDict) 

481 self.log.info("Selected PSF sizes:\n%s", psfSizeDict) 

482 self.bufferSize = int(np.ceil(np.max(dcrShifts)) + 1) 

483 # Adding dictionary directly to the metadata did not add all elements. 

484 # Ensure that the metadata is listed in the correct order. 

485 visits = airmassDict.keys() 

486 self.metadata["visits"] = list(visits) 

487 self.metadata["airmasses"] = [airmassDict[v] for v in visits] 

488 self.metadata["parallacticAngles"] = [angleDict[v] for v in visits] 

489 self.metadata["selectedPsfSizes"] = [psfSizeDict[v] for v in visits] 

490 try: 

491 psf = self.selectCoaddPsf(templateCoadd, warpRefList) 

492 except Exception as e: 

493 self.log.warning("Unable to calculate restricted PSF, using default coadd PSF: %s", e) 

494 else: 

495 psf = templateCoadd.getPsf() 

496 dcrModels = DcrModel.fromImage( 

497 templateCoadd.maskedImage, 

498 self.config.dcrNumSubfilters, 

499 effectiveWavelength=self.config.effectiveWavelength, 

500 bandwidth=self.config.bandwidth, 

501 wcs=templateCoadd.getWcs(), 

502 filterLabel=filterLabel, 

503 psf=psf, 

504 ) 

505 return dcrModels 

506 

507 @timeMethod 

508 def run(self, skyInfo, *, warpRefList, imageScalerList, weightList, supplementaryData=None, **kwargs): 

509 r"""Assemble the coadd. 

510 

511 Requires additional inputs Struct ``supplementaryData`` to contain a 

512 ``templateCoadd`` that serves as the model of the static sky. 

513 

514 Find artifacts and apply them to the warps' masks creating a list of 

515 alternative masks with a new "CLIPPED" plane and updated "NO_DATA" 

516 plane then pass these alternative masks to the base class's assemble 

517 method. 

518 

519 Divide the ``templateCoadd`` evenly between each subfilter of a 

520 ``DcrModel`` as the starting best estimate of the true wavelength- 

521 dependent sky. Forward model the ``DcrModel`` using the known 

522 chromatic effects in each subfilter and calculate a convergence metric 

523 based on how well the modeled template matches the input warps. If 

524 the convergence has not yet reached the desired threshold, then shift 

525 and stack the residual images to build a new ``DcrModel``. Apply 

526 conditioning to prevent oscillating solutions between iterations or 

527 between subfilters. 

528 

529 Once the ``DcrModel`` reaches convergence or the maximum number of 

530 iterations has been reached, fill the metadata for each subfilter 

531 image and make them proper ``coaddExposure``\ s. 

532 

533 Parameters 

534 ---------- 

535 skyInfo : `lsst.pipe.base.Struct` 

536 Patch geometry information, from getSkyInfo 

537 warpRefList : `list` [`lsst.daf.butler.DeferredDatasetHandle`] 

538 The data references to the input warped exposures. 

539 imageScalerList : `list` [`lsst.pipe.task.ImageScaler`] 

540 The image scalars correct for the zero point of the exposures. 

541 Deprecated and will be removed after v29 in DM-49083. 

542 weightList : `list` [`float`] 

543 The weight to give each input exposure in the coadd. 

544 supplementaryData : `lsst.pipe.base.Struct` 

545 Result struct returned by ``_makeSupplementaryData`` with 

546 attributes: 

547 

548 ``templateCoadd`` 

549 Coadded exposure (`lsst.afw.image.Exposure`). 

550 

551 Returns 

552 ------- 

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

554 Results as a struct with attributes: 

555 

556 ``coaddExposure`` 

557 Coadded exposure (`lsst.afw.image.Exposure`). 

558 ``nImage`` 

559 Exposure count image (`lsst.afw.image.ImageU`). 

560 ``dcrCoadds`` 

561 `list` of coadded exposures for each subfilter. 

562 ``dcrNImages`` 

563 `list` of exposure count images for each subfilter. 

564 """ 

565 minNumIter = self.config.minNumIter or self.config.dcrNumSubfilters 

566 maxNumIter = self.config.maxNumIter or self.config.dcrNumSubfilters * 3 

567 templateCoadd = supplementaryData.templateCoadd 

568 baseMask = templateCoadd.mask.clone() 

569 # The variance plane is for each subfilter 

570 # and should be proportionately lower than the full-band image 

571 baseVariance = templateCoadd.variance.clone() 

572 baseVariance /= self.config.dcrNumSubfilters 

573 spanSetMaskList = self.findArtifacts(templateCoadd, warpRefList, imageScalerList) 

574 # Note that the mask gets cleared in ``findArtifacts``, but we want to 

575 # preserve the mask. 

576 templateCoadd.setMask(baseMask) 

577 badMaskPlanes = self.config.badMaskPlanes[:] 

578 # Note that is important that we do not add "CLIPPED" to 

579 # ``badMaskPlanes``. This is because pixels in observations that are 

580 # significantly affected by DCR are likely to have many pixels that are 

581 # both "DETECTED" and "CLIPPED", but those are necessary to constrain 

582 # the DCR model. 

583 badPixelMask = templateCoadd.mask.getPlaneBitMask(badMaskPlanes) 

584 

585 stats = self.prepareStats(mask=badPixelMask) 

586 dcrModels = self.prepareDcrInputs(templateCoadd, warpRefList, weightList) 

587 if self.config.doNImage: 

588 dcrNImages, dcrWeights = self.calculateNImage( 

589 dcrModels, skyInfo.bbox, warpRefList, spanSetMaskList, stats.ctrl 

590 ) 

591 nImage = afwImage.ImageU(skyInfo.bbox) 

592 # Note that this nImage will be a factor of dcrNumSubfilters higher 

593 # than the nImage returned by assembleCoadd for most pixels. This 

594 # is because each subfilter may have a different nImage, and 

595 # fractional values are not allowed. 

596 for dcrNImage in dcrNImages: 

597 nImage += dcrNImage 

598 else: 

599 dcrNImages = None 

600 

601 subregionSize = geom.Extent2I(*self.config.subregionSize) 

602 nSubregions = ceil(skyInfo.bbox.getHeight() / subregionSize[1]) * ceil( 

603 skyInfo.bbox.getWidth() / subregionSize[0] 

604 ) 

605 subIter = 0 

606 initialConvergence = self.calculateConvergence( 

607 dcrModels, None, skyInfo.bbox, warpRefList, weightList, stats.ctrl 

608 ) 

609 self.metadata["initialConvergence"] = initialConvergence 

610 self.log.info("Initial convergence for full patch %.4f%%", initialConvergence) 

611 for subBBox in subBBoxIter(skyInfo.bbox, subregionSize): 

612 modelIter = 0 

613 subIter += 1 

614 self.log.info( 

615 "Computing coadd over patch %s subregion %s of %s: %s", 

616 skyInfo.patchInfo.getIndex(), 

617 subIter, 

618 nSubregions, 

619 subBBox, 

620 ) 

621 dcrBBox = geom.Box2I(subBBox) 

622 dcrBBox.grow(self.bufferSize) 

623 dcrBBox.clip(dcrModels.bbox) 

624 modelWeights = self.calculateModelWeights(dcrModels, dcrBBox) 

625 subExposures = self.loadSubExposures( 

626 dcrBBox, stats.ctrl, warpRefList, imageScalerList, spanSetMaskList 

627 ) 

628 convergenceMetric = self.calculateConvergence( 

629 dcrModels, subExposures, subBBox, warpRefList, weightList, stats.ctrl 

630 ) 

631 self.log.info("Initial convergence : %s", convergenceMetric) 

632 convergenceList = [convergenceMetric] 

633 gainList = [] 

634 convergenceCheck = 1.0 

635 refImage = templateCoadd.image 

636 while convergenceCheck > self.config.convergenceThreshold or modelIter <= minNumIter: 

637 gain = self.calculateGain(convergenceList, gainList) 

638 self.dcrAssembleSubregion( 

639 dcrModels, 

640 subExposures, 

641 subBBox, 

642 dcrBBox, 

643 warpRefList, 

644 stats.ctrl, 

645 convergenceMetric, 

646 gain, 

647 modelWeights, 

648 refImage, 

649 dcrWeights, 

650 ) 

651 if self.config.useConvergence: 

652 convergenceMetric = self.calculateConvergence( 

653 dcrModels, subExposures, subBBox, warpRefList, weightList, stats.ctrl 

654 ) 

655 if convergenceMetric == 0: 

656 self.log.warning( 

657 "Coadd patch %s subregion %s had convergence metric of 0.0 which is " 

658 "most likely due to there being no valid data in the region.", 

659 skyInfo.patchInfo.getIndex(), 

660 subIter, 

661 ) 

662 break 

663 convergenceCheck = (convergenceList[-1] - convergenceMetric) / convergenceMetric 

664 if (convergenceCheck < 0) & (modelIter > minNumIter): 

665 self.log.warning( 

666 "Coadd patch %s subregion %s diverged before reaching maximum " 

667 "iterations or desired convergence improvement of %s." 

668 " Divergence: %s", 

669 skyInfo.patchInfo.getIndex(), 

670 subIter, 

671 self.config.convergenceThreshold, 

672 convergenceCheck, 

673 ) 

674 break 

675 convergenceList.append(convergenceMetric) 

676 if modelIter > maxNumIter: 

677 if self.config.useConvergence: 

678 self.log.warning( 

679 "Coadd patch %s subregion %s reached maximum iterations " 

680 "before reaching desired convergence improvement of %s." 

681 " Final convergence improvement: %s", 

682 skyInfo.patchInfo.getIndex(), 

683 subIter, 

684 self.config.convergenceThreshold, 

685 convergenceCheck, 

686 ) 

687 break 

688 

689 if self.config.useConvergence: 

690 self.log.info( 

691 "Iteration %s with convergence metric %s, %.4f%% improvement (gain: %.2f)", 

692 modelIter, 

693 convergenceMetric, 

694 100.0 * convergenceCheck, 

695 gain, 

696 ) 

697 modelIter += 1 

698 else: 

699 if self.config.useConvergence: 

700 self.log.info( 

701 "Coadd patch %s subregion %s finished with " 

702 "convergence metric %s after %s iterations", 

703 skyInfo.patchInfo.getIndex(), 

704 subIter, 

705 convergenceMetric, 

706 modelIter, 

707 ) 

708 else: 

709 self.log.info( 

710 "Coadd patch %s subregion %s finished after %s iterations", 

711 skyInfo.patchInfo.getIndex(), 

712 subIter, 

713 modelIter, 

714 ) 

715 if self.config.useConvergence and convergenceMetric > 0: 

716 self.log.info( 

717 "Final convergence improvement was %.4f%% overall", 

718 100 * (convergenceList[0] - convergenceMetric) / convergenceMetric, 

719 ) 

720 

721 finalConvergence = self.calculateConvergence( 

722 dcrModels, None, skyInfo.bbox, warpRefList, weightList, stats.ctrl 

723 ) 

724 self.metadata["finalConvergence"] = finalConvergence 

725 self.log.info("Final convergence for full patch %.4f%%", finalConvergence) 

726 

727 # Remove in DM-49083 

728 if self.config.doScaleZeroPoint: 

729 calibration = self.scaleZeroPoint.getPhotoCalib() 

730 else: 

731 calibration = afwImage.PhotoCalib(1.0) 

732 

733 dcrCoadds = self.fillCoadd( 

734 dcrModels, 

735 skyInfo, 

736 warpRefList, 

737 weightList, 

738 calibration=calibration, 

739 coaddInputs=templateCoadd.getInfo().getCoaddInputs(), 

740 mask=baseMask, 

741 variance=baseVariance, 

742 ) 

743 coaddExposure = self.stackCoadd(dcrCoadds) 

744 return pipeBase.Struct( 

745 coaddExposure=coaddExposure, nImage=nImage, dcrCoadds=dcrCoadds, dcrNImages=dcrNImages 

746 ) 

747 

748 def calculateNImage(self, dcrModels, bbox, warpRefList, spanSetMaskList, statsCtrl): 

749 """Calculate the number of exposures contributing to each subfilter. 

750 

751 Parameters 

752 ---------- 

753 dcrModels : `lsst.pipe.tasks.DcrModel` 

754 Best fit model of the true sky after correcting chromatic effects. 

755 bbox : `lsst.geom.box.Box2I` 

756 Bounding box of the patch to coadd. 

757 warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` 

758 The data references to the input warped exposures. 

759 spanSetMaskList : `list` of `dict` containing spanSet lists, or `None` 

760 Each element of the `dict` contains the new mask plane name 

761 (e.g. "CLIPPED and/or "NO_DATA") as the key, 

762 and the list of SpanSets to apply to the mask. 

763 statsCtrl : `lsst.afw.math.StatisticsControl` 

764 Statistics control object for coadd 

765 

766 Returns 

767 ------- 

768 dcrNImages : `list` of `lsst.afw.image.ImageU` 

769 List of exposure count images for each subfilter. 

770 dcrWeights : `list` of `lsst.afw.image.ImageF` 

771 Per-pixel weights for each subfilter. 

772 Equal to 1/(number of unmasked images contributing to each pixel). 

773 """ 

774 dcrNImages = [afwImage.ImageU(bbox) for subfilter in range(self.config.dcrNumSubfilters)] 

775 dcrWeights = [afwImage.ImageF(bbox) for subfilter in range(self.config.dcrNumSubfilters)] 

776 for warpExpRef, altMaskSpans in zip(warpRefList, spanSetMaskList): 

777 exposure = warpExpRef.get(parameters={"bbox": bbox}) 

778 visitInfo = exposure.getInfo().getVisitInfo() 

779 wcs = exposure.getInfo().getWcs() 

780 mask = exposure.mask 

781 if altMaskSpans is not None: 

782 self.applyAltMaskPlanes(mask, altMaskSpans) 

783 weightImage = np.zeros_like(exposure.image.array) 

784 weightImage[(mask.array & statsCtrl.getAndMask()) == 0] = 1.0 

785 # The weights must be shifted in exactly the same way as the 

786 # residuals, because they will be used as the denominator in the 

787 # weighted average of residuals. 

788 weightsGenerator = self.dcrResiduals( 

789 weightImage, visitInfo, wcs, dcrModels.effectiveWavelength, dcrModels.bandwidth 

790 ) 

791 for shiftedWeights, dcrNImage, dcrWeight in zip(weightsGenerator, dcrNImages, dcrWeights): 

792 dcrNImage.array += np.rint(shiftedWeights).astype(dcrNImage.array.dtype) 

793 dcrWeight.array += shiftedWeights 

794 # Exclude any pixels that don't have at least one exposure contributing 

795 # in all subfilters 

796 weightsThreshold = 1.0 

797 goodPix = dcrWeights[0].array > weightsThreshold 

798 for weights in dcrWeights[1:]: 

799 goodPix = (weights.array > weightsThreshold) & goodPix 

800 for subfilter in range(self.config.dcrNumSubfilters): 

801 dcrWeights[subfilter].array[goodPix] = 1.0 / dcrWeights[subfilter].array[goodPix] 

802 dcrWeights[subfilter].array[~goodPix] = 0.0 

803 dcrNImages[subfilter].array[~goodPix] = 0 

804 return (dcrNImages, dcrWeights) 

805 

806 def dcrAssembleSubregion( 

807 self, 

808 dcrModels, 

809 subExposures, 

810 bbox, 

811 dcrBBox, 

812 warpRefList, 

813 statsCtrl, 

814 convergenceMetric, 

815 gain, 

816 modelWeights, 

817 refImage, 

818 dcrWeights, 

819 ): 

820 """Assemble the DCR coadd for a sub-region. 

821 

822 Build a DCR-matched template for each input exposure, then shift the 

823 residuals according to the DCR in each subfilter. 

824 Stack the shifted residuals and apply them as a correction to the 

825 solution from the previous iteration. 

826 Restrict the new model solutions from varying by more than a factor of 

827 `modelClampFactor` from the last solution, and additionally restrict 

828 the individual subfilter models from varying by more than a factor of 

829 `frequencyClampFactor` from their average. 

830 Finally, mitigate potentially oscillating solutions by averaging the 

831 new solution with the solution from the previous iteration, weighted by 

832 their convergence metric. 

833 

834 Parameters 

835 ---------- 

836 dcrModels : `lsst.pipe.tasks.DcrModel` 

837 Best fit model of the true sky after correcting chromatic effects. 

838 subExposures : `dict` of `lsst.afw.image.ExposureF` 

839 The pre-loaded exposures for the current subregion. 

840 bbox : `lsst.geom.box.Box2I` 

841 Bounding box of the subregion to coadd. 

842 dcrBBox : `lsst.geom.box.Box2I` 

843 Sub-region of the coadd which includes a buffer to allow for DCR. 

844 warpRefList : `list` of `lsst.daf.butler.DeferredDatasetHandle` 

845 The data references to the input warped exposures. 

846 statsCtrl : `lsst.afw.math.StatisticsControl` 

847 Statistics control object for coadd. 

848 convergenceMetric : `float` 

849 Quality of fit metric for the matched templates of the input 

850 images. 

851 gain : `float`, optional 

852 Relative weight to give the new solution when updating the model. 

853 modelWeights : `numpy.ndarray` or `float` 

854 A 2D array of weight values that tapers smoothly to zero away from 

855 detected sources. Set to a placeholder value of 1.0 if 

856 ``self.config.useModelWeights`` is False. 

857 refImage : `lsst.afw.image.Image` 

858 A reference image used to supply the default pixel values. 

859 dcrWeights : `list` of `lsst.afw.image.Image` 

860 Per-pixel weights for each subfilter. 

861 Equal to 1/(number of unmasked images contributing to each pixel). 

862 """ 

863 residualGeneratorList = [] 

864 

865 for warpExpRef in warpRefList: 

866 visit = warpExpRef.dataId["visit"] 

867 exposure = subExposures[visit] 

868 visitInfo = exposure.getInfo().getVisitInfo() 

869 wcs = exposure.getInfo().getWcs() 

870 templateImage = dcrModels.buildMatchedTemplate( 

871 exposure=exposure, 

872 bbox=exposure.getBBox(), 

873 order=self.config.imageInterpOrder, 

874 splitSubfilters=self.config.splitSubfilters, 

875 splitThreshold=self.config.splitThreshold, 

876 amplifyModel=self.config.accelerateModel, 

877 ) 

878 residual = exposure.image.array - templateImage.array 

879 # Note that the variance plane here is used to store weights, not 

880 # the actual variance 

881 residual *= exposure.variance.array 

882 # The residuals are stored as a list of generators. 

883 # This allows the residual for a given subfilter and exposure to be 

884 # created on the fly, instead of needing to store them all in 

885 # memory. 

886 residualGeneratorList.append( 

887 self.dcrResiduals( 

888 residual, visitInfo, wcs, dcrModels.effectiveWavelength, dcrModels.bandwidth 

889 ) 

890 ) 

891