Coverage for python / lsst / drp / tasks / make_direct_warp.py: 24%

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 

22from __future__ import annotations 

23 

24import dataclasses 

25from collections.abc import Mapping 

26from typing import TYPE_CHECKING, Iterable 

27 

28import numpy as np 

29 

30from lsst.afw.image import ExposureF, Image, Mask, PhotoCalib 

31from lsst.afw.math import BackgroundList, Warper 

32from lsst.coadd.utils import copyGoodPixels 

33from lsst.daf.butler import DataCoordinate, DeferredDatasetHandle 

34from lsst.geom import Box2D 

35from lsst.meas.algorithms import CoaddPsf, CoaddPsfConfig, backgroundFlatContext 

36from lsst.meas.algorithms.cloughTocher2DInterpolator import CloughTocher2DInterpolateTask 

37from lsst.meas.base import DetectorVisitIdGeneratorConfig 

38from lsst.pex.config import ConfigField, ConfigurableField, Field, RangeField 

39from lsst.pipe.base import ( 

40 InMemoryDatasetHandle, 

41 NoWorkFound, 

42 PipelineTask, 

43 PipelineTaskConfig, 

44 PipelineTaskConnections, 

45 Struct, 

46) 

47from lsst.pipe.base.connectionTypes import Input, Output 

48from lsst.pipe.tasks.coaddBase import makeSkyInfo 

49from lsst.pipe.tasks.coaddInputRecorder import CoaddInputRecorderTask 

50from lsst.pipe.tasks.selectImages import PsfWcsSelectImagesTask 

51from lsst.skymap import BaseSkyMap 

52 

53if TYPE_CHECKING: 

54 from lsst.afw.image import MaskedImage 

55 from lsst.afw.table import ExposureCatalog 

56 

57 

58__all__ = ( 

59 "MakeDirectWarpConfig", 

60 "MakeDirectWarpTask", 

61) 

62 

63 

64@dataclasses.dataclass 

65class WarpDetectorInputs: 

66 """Inputs passed to `MakeDirectWarpTask.run` for a single detector.""" 

67 

68 exposure_or_handle: ExposureF | DeferredDatasetHandle | InMemoryDatasetHandle 

69 """Detector image with initial calibration objects, or a deferred-load 

70 handle for one. 

71 """ 

72 

73 data_id: DataCoordinate 

74 """Butler data ID for this detector.""" 

75 

76 background_revert: BackgroundList | None = None 

77 """Background model to restore in (i.e. add to) the image.""" 

78 

79 background_apply: BackgroundList | None = None 

80 """Background model to apply to (i.e. subtract from) the image.""" 

81 

82 background_ratio_or_handle: Image | DeferredDatasetHandle | InMemoryDatasetHandle | None = None 

83 """Ratio of background-flattened image to photometric-flattened image.""" 

84 

85 @property 

86 def exposure(self) -> ExposureF: 

87 """Get the exposure object, loading it if necessary.""" 

88 if not isinstance(self.exposure_or_handle, ExposureF): 

89 self.exposure_or_handle = self.exposure_or_handle.get() 

90 

91 return self.exposure_or_handle 

92 

93 @property 

94 def background_to_photometric_ratio(self) -> Image: 

95 """Get the background_to_photometric object, loading if necessary.""" 

96 if isinstance(self.background_ratio_or_handle, (DeferredDatasetHandle, InMemoryDatasetHandle)): 

97 self.background_ratio_or_handle = self.background_ratio_or_handle.get() 

98 

99 return self.background_ratio_or_handle 

100 

101 def apply_background(self) -> None: 

102 """Apply (subtract) the `background_apply` to the exposure in-place. 

103 

104 Raises 

105 ------ 

106 RuntimeError 

107 Raised if `background_apply` is None. 

108 """ 

109 if self.background_apply is None: 

110 raise RuntimeError("No background to apply") 

111 

112 if self.background_apply: 

113 with backgroundFlatContext( 

114 self.exposure.maskedImage, 

115 self.background_to_photometric_ratio is not None, 

116 backgroundToPhotometricRatio=self.background_to_photometric_ratio, 

117 ): 

118 self.exposure.maskedImage -= self.background_apply.getImage() 

119 

120 def revert_background(self) -> None: 

121 """Revert (add) the `background_revert` from the exposure in-place. 

122 

123 Raises 

124 ------ 

125 RuntimeError 

126 Raised if `background_revert` is None. 

127 """ 

128 if self.background_revert is None: 

129 raise RuntimeError("No background to revert") 

130 

131 if self.background_revert: 

132 with backgroundFlatContext( 

133 self.exposure.maskedImage, 

134 self.background_to_photometric_ratio is not None, 

135 backgroundToPhotometricRatio=self.background_to_photometric_ratio, 

136 ): 

137 # Add only if `background_revert` is not a trivial background. 

138 self.exposure.maskedImage += self.background_revert.getImage() 

139 

140 

141class MakeDirectWarpConnections( 

142 PipelineTaskConnections, 

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

144 defaultTemplates={ 

145 "coaddName": "deep", 

146 }, 

147): 

148 """Connections for MakeWarpTask""" 

149 

150 calexp_list = Input( 

151 doc="Input exposures to be interpolated and resampled onto a SkyMap " "projection/patch.", 

152 name="calexp", 

153 storageClass="ExposureF", 

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

155 multiple=True, 

156 deferLoad=True, 

157 ) 

158 background_revert_list = Input( 

159 doc="Background to be reverted (i.e., added back to the calexp). " 

160 "This connection is used only if doRevertOldBackground=True.", 

161 name="calexpBackground", 

162 storageClass="Background", 

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

164 multiple=True, 

165 ) 

166 background_apply_list = Input( 

167 doc="Background to be applied (subtracted from the calexp). " 

168 "This is used only if doApplyNewBackground=True.", 

169 name="skyCorr", 

170 storageClass="Background", 

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

172 multiple=True, 

173 ) 

174 background_to_photometric_ratio_list = Input( 

175 doc="Ratio of a background-flattened image to a photometric-flattened image. " 

176 "This is used only if doRevertOldBackground=True or doApplyNewBackground=True " 

177 "and doApplyFlatBackgroundRatio=True.", 

178 name="background_to_photometric_ratio", 

179 storageClass="Image", 

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

181 multiple=True, 

182 deferLoad=True, 

183 ) 

184 visit_summary = Input( 

185 doc="Input visit-summary catalog with updated calibration objects.", 

186 name="finalVisitSummary", 

187 storageClass="ExposureCatalog", 

188 dimensions=("instrument", "visit"), 

189 ) 

190 sky_map = Input( 

191 doc="Input definition of geometry/bbox and projection/wcs for warps.", 

192 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

193 storageClass="SkyMap", 

194 dimensions=("skymap",), 

195 ) 

196 # Declare all possible outputs (except noise, which is configurable) 

197 warp = Output( 

198 doc="Output direct warped exposure produced by resampling calexps " "onto the skyMap patch geometry.", 

199 name="{coaddName}Coadd_directWarp", 

200 storageClass="ExposureF", 

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

202 ) 

203 masked_fraction_warp = Output( 

204 doc="Output masked fraction warped exposure.", 

205 name="{coaddName}Coadd_directWarp_maskedFraction", 

206 storageClass="ImageF", 

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

208 ) 

209 

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

211 super().__init__(config=config) 

212 if not config: 

213 return 

214 

215 if not config.doRevertOldBackground: 

216 del self.background_revert_list 

217 if not config.doApplyNewBackground: 

218 del self.background_apply_list 

219 if not config.doApplyFlatBackgroundRatio or ( 

220 not config.doRevertOldBackground and not config.doApplyNewBackground 

221 ): 

222 del self.background_to_photometric_ratio_list 

223 

224 if not config.doWarpMaskedFraction: 

225 del self.masked_fraction_warp 

226 

227 # Dynamically set output connections for noise images, depending on the 

228 # number of noise realization specified in the config. 

229 for n in range(config.numberOfNoiseRealizations): 

230 noise_warp = Output( 

231 doc=f"Output direct warped noise exposure ({n})", 

232 name=f"direct_warp_noise{n}", 

233 # Store it as a MaskedImage to preserve the variance plane. 

234 storageClass="MaskedImageF", 

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

236 ) 

237 setattr(self, f"noise_warp{n}", noise_warp) 

238 

239 

240class MakeDirectWarpConfig( 

241 PipelineTaskConfig, 

242 pipelineConnections=MakeDirectWarpConnections, 

243): 

244 """Configuration for the MakeDirectWarpTask. 

245 

246 The config fields are as similar as possible to the corresponding fields in 

247 MakeWarpConfig. 

248 

249 Notes 

250 ----- 

251 The config fields are in camelCase to match the fields in the earlier 

252 version of the makeWarp task as closely as possible. 

253 """ 

254 

255 MAX_NUMBER_OF_NOISE_REALIZATIONS = 3 

256 """ 

257 numberOfNoiseRealizations is defined as a RangeField to prevent from 

258 making multiple output connections and blowing up the memory usage by 

259 accident. An upper bound of 3 is based on the best guess of the maximum 

260 number of noise realizations that will be used for metadetection. 

261 """ 

262 

263 numberOfNoiseRealizations = RangeField[int]( 

264 doc="Number of noise realizations to simulate and persist.", 

265 default=0, 

266 min=0, 

267 max=MAX_NUMBER_OF_NOISE_REALIZATIONS, 

268 inclusiveMax=True, 

269 ) 

270 seedOffset = Field[int]( 

271 doc="Offset to the seed used for the noise realization. This can be " 

272 "used to create a different noise realization if the default ones " 

273 "are catastrophic, or for testing sensitivity to the noise.", 

274 default=0, 

275 ) 

276 useMedianVariance = Field[bool]( 

277 doc="Use the median of variance plane in the input calexp to generate " 

278 "noise realizations? If False, per-pixel variance will be used.", 

279 default=True, 

280 ) 

281 doRevertOldBackground = Field[bool]( 

282 doc="Revert the old backgrounds from the `background_revert_list` " "connection?", 

283 default=False, 

284 ) 

285 doApplyNewBackground = Field[bool]( 

286 doc="Apply the new backgrounds from the `background_apply_list` " "connection?", 

287 default=False, 

288 ) 

289 doApplyFlatBackgroundRatio = Field[bool]( 

290 doc="Apply flat background ratio prior to background adjustments? Should be True " 

291 "if processing was done with an illumination correction.", 

292 default=False, 

293 ) 

294 useVisitSummaryPsf = Field[bool]( 

295 doc="If True, use the PSF model and aperture corrections from the " 

296 "'visit_summary' connection to make the warp. If False, use the " 

297 "PSF model and aperture corrections from the 'calexp' connection.", 

298 default=True, 

299 ) 

300 useVisitSummaryWcs = Field[bool]( 

301 doc="If True, use the WCS from the " 

302 "'visit_summary' connection to make the warp. If False, use the " 

303 "WCS from the 'calexp' connection.", 

304 default=True, 

305 ) 

306 useVisitSummaryPhotoCalib = Field[bool]( 

307 doc="If True, use the photometric calibration from the " 

308 "'visit_summary' connection to make the warp. If False, use the " 

309 "photometric calibration from the 'calexp' connection.", 

310 default=True, 

311 ) 

312 doSelectPreWarp = Field[bool]( 

313 doc="Select ccds before warping?", 

314 default=True, 

315 ) 

316 select = ConfigurableField( 

317 doc="Image selection subtask.", 

318 target=PsfWcsSelectImagesTask, 

319 ) 

320 doPreWarpInterpolation = Field[bool]( 

321 doc="Interpolate over bad pixels before warping?", 

322 default=False, 

323 ) 

324 preWarpInterpolation = ConfigurableField( 

325 doc="Interpolation task to use for pre-warping interpolation", 

326 target=CloughTocher2DInterpolateTask, 

327 ) 

328 inputRecorder = ConfigurableField( 

329 doc="Subtask that helps fill CoaddInputs catalogs added to the final " "coadd", 

330 target=CoaddInputRecorderTask, 

331 ) 

332 includeCalibVar = Field[bool]( 

333 doc="Add photometric calibration variance to warp variance plane?", 

334 default=False, 

335 deprecated="Deprecated and disabled. Will be removed after v29.", 

336 ) 

337 border = Field[int]( 

338 doc="Pad the patch boundary of the warp by these many pixels, so as to allow for PSF-matching later", 

339 default=256, 

340 ) 

341 warper = ConfigField( 

342 doc="Configuration for the warper that warps the image and noise", 

343 dtype=Warper.ConfigClass, 

344 ) 

345 doWarpMaskedFraction = Field[bool]( 

346 doc="Warp the masked fraction image?", 

347 default=True, 

348 ) 

349 maskedFractionWarper = ConfigField( 

350 doc="Configuration for the warp that warps the mask fraction image", 

351 dtype=Warper.ConfigClass, 

352 ) 

353 coaddPsf = ConfigField( 

354 doc="Configuration for CoaddPsf", 

355 dtype=CoaddPsfConfig, 

356 ) 

357 idGenerator = DetectorVisitIdGeneratorConfig.make_field() 

358 

359 # Use bgSubtracted and doApplySkyCorr to match the old MakeWarpConfig, 

360 # but as properties instead of config fields. 

361 @property 

362 def bgSubtracted(self) -> bool: 

363 return not self.doRevertOldBackground 

364 

365 @bgSubtracted.setter 

366 def bgSubtracted(self, value: bool) -> None: 

367 self.doRevertOldBackground = not value 

368 

369 @property 

370 def doApplySkyCorr(self) -> bool: 

371 return self.doApplyNewBackground 

372 

373 @doApplySkyCorr.setter 

374 def doApplySkyCorr(self, value: bool) -> None: 

375 self.doApplyNewBackground = value 

376 

377 def setDefaults(self) -> None: 

378 super().setDefaults() 

379 self.warper.warpingKernelName = "lanczos3" 

380 self.warper.cacheSize = 0 

381 self.maskedFractionWarper.warpingKernelName = "bilinear" 

382 

383 

384class MakeDirectWarpTask(PipelineTask): 

385 """Warp single-detector images onto a common projection. 

386 

387 This task iterates over multiple images (corresponding to different 

388 detectors) from a single visit that overlap the target patch. Pixels that 

389 receive no input from any detector are set to NaN in the output image, and 

390 NO_DATA bit is set in the mask plane. 

391 

392 This differs from the standard `MakeWarp` Task in the following 

393 ways: 

394 

395 1. No selection on ccds at the time of warping. This is done later during 

396 the coaddition stage. 

397 2. Interpolate over a set of masked pixels before warping. 

398 3. Generate an image where each pixel denotes how much of the pixel is 

399 masked. 

400 4. Generate multiple noise warps with the same interpolation applied. 

401 5. No option to produce a PSF-matched warp. 

402 """ 

403 

404 ConfigClass = MakeDirectWarpConfig 

405 _DefaultName = "makeDirectWarp" 

406 

407 def __init__(self, **kwargs): 

408 super().__init__(**kwargs) 

409 self.makeSubtask("inputRecorder") 

410 self.makeSubtask("preWarpInterpolation") 

411 if self.config.doSelectPreWarp: 

412 self.makeSubtask("select") 

413 

414 self.warper = Warper.fromConfig(self.config.warper) 

415 if self.config.doWarpMaskedFraction: 

416 self.maskedFractionWarper = Warper.fromConfig(self.config.maskedFractionWarper) 

417 

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

419 # Docstring inherited. 

420 

421 inputs: Mapping[int, WarpDetectorInputs] = {} # Detector ID -> WarpDetectorInputs 

422 for handle in butlerQC.get(inputRefs.calexp_list): 

423 inputs[handle.dataId["detector"]] = WarpDetectorInputs( 

424 exposure_or_handle=handle, 

425 data_id=handle.dataId, 

426 ) 

427 

428 if not inputs: 

429 raise NoWorkFound("No input warps provided for co-addition") 

430 

431 # Add backgrounds to the inputs struct, being careful not to assume 

432 # they're present for the same detectors we got image handles for, in 

433 # case of upstream errors. 

434 for ref in getattr(inputRefs, "background_revert_list", []): 

435 inputs[ref.dataId["detector"]].background_revert = butlerQC.get(ref) 

436 for ref in getattr(inputRefs, "background_apply_list", []): 

437 inputs[ref.dataId["detector"]].background_apply = butlerQC.get(ref) 

438 for ref in getattr(inputRefs, "background_to_photometric_ratio_list", []): 

439 inputs[ref.dataId["detector"]].background_ratio_or_handle = butlerQC.get(ref) 

440 

441 visit_summary = butlerQC.get(inputRefs.visit_summary) 

442 sky_map = butlerQC.get(inputRefs.sky_map) 

443 

444 quantumDataId = butlerQC.quantum.dataId 

445 sky_info = makeSkyInfo( 

446 sky_map, 

447 tractId=quantumDataId["tract"], 

448 patchId=quantumDataId["patch"], 

449 ) 

450 

451 results = self.run(inputs, sky_info, visit_summary=visit_summary) 

452 butlerQC.put(results, outputRefs) 

453 

454 def _preselect_inputs( 

455 self, 

456 inputs: Mapping[int, WarpDetectorInputs], 

457 sky_info: Struct, 

458 visit_summary: ExposureCatalog, 

459 ) -> dict[int, WarpDetectorInputs]: 

460 """Filter the list of inputs via a 'select' subtask. 

461 

462 Parameters 

463 ---------- 

464 inputs : ``~collections.abc.Mapping` [ `int`, `WarpDetectorInputs` ] 

465 Per-detector input structs. 

466 sky_info : `lsst.pipe.base.Struct` 

467 Structure with information about the tract and patch. 

468 visit_summary : `~lsst.afw.table.ExposureCatalog` 

469 Record with updated calibration information for the full visit. 

470 

471 Returns 

472 ------- 

473 filtered_inputs : `dict` [ `int`, `WarpDetectorInputs` ] 

474 Like ``inputs``, with rejected detectors dropped. 

475 """ 

476 data_id_list, bbox_list, wcs_list = [], [], [] 

477 for detector_id, detector_inputs in inputs.items(): 

478 row = visit_summary.find(detector_id) 

479 if row is None: 

480 self.log.info( 

481 "Input calexp is not listed in visit_summary: %d; assuming bad.", 

482 detector_id, 

483 ) 

484 continue 

485 data_id_list.append(detector_inputs.data_id) 

486 bbox_list.append(row.getBBox()) 

487 wcs_list.append(row.getWcs()) 

488 

489 cornerPosList = Box2D(sky_info.bbox).getCorners() 

490 coordList = [sky_info.wcs.pixelToSky(pos) for pos in cornerPosList] 

491 

492 good_indices = self.select.run( 

493 bboxList=bbox_list, 

494 wcsList=wcs_list, 

495 visitSummary=visit_summary, 

496 coordList=coordList, 

497 dataIds=data_id_list, 

498 ) 

499 detector_ids = list(inputs.keys()) 

500 good_detector_ids = [detector_ids[i] for i in good_indices] 

501 return {detector_id: inputs[detector_id] for detector_id in good_detector_ids} 

502 

503 def run(self, inputs: Mapping[int, WarpDetectorInputs], sky_info, visit_summary) -> Struct: 

504 """Create a Warp dataset from inputs. 

505 

506 Parameters 

507 ---------- 

508 inputs : `Mapping` [ `int`, `WarpDetectorInputs` ] 

509 Dictionary of input datasets, with the detector id being the key. 

510 sky_info : `~lsst.pipe.base.Struct` 

511 A Struct object containing wcs, bounding box, and other information 

512 about the patches within the tract. 

513 visit_summary : `~lsst.afw.table.ExposureCatalog` | None 

514 Table of visit summary information. If provided, the visit summary 

515 information will be used to update the calibration of the input 

516 exposures. If None, the input exposures will be used as-is. 

517 

518 Returns 

519 ------- 

520 results : `~lsst.pipe.base.Struct` 

521 A Struct object containing the warped exposure, noise exposure(s), 

522 and masked fraction image. 

523 """ 

524 if self.config.doSelectPreWarp: 

525 inputs = self._preselect_inputs(inputs, sky_info, visit_summary=visit_summary) 

526 if not inputs: 

527 raise NoWorkFound("No input warps remain after selection for co-addition") 

528 

529 sky_info.bbox.grow(self.config.border) 

530 target_bbox, target_wcs = sky_info.bbox, sky_info.wcs 

531 

532 # Initialize the objects that will hold the warp. 

533 final_warp = ExposureF(target_bbox, target_wcs) 

534 

535 for _, warp_detector_input in inputs.items(): 

536 visit_id = warp_detector_input.data_id["visit"] 

537 break # Just need the visit id from any one of the inputs. 

538 

539 # The warpExposure routine is expensive, and we do not want to call 

540 # it twice (i.e., a second time for PSF-matched warps). We do not 

541 # want to hold all the warped exposures in memory at once either. 

542 # So we create empty exposure(s) to accumulate the warps of each type, 

543 # and then process each detector serially. 

544 final_warp = self._prepareEmptyExposure(sky_info) 

545 final_masked_fraction_warp = self._prepareEmptyExposure(sky_info) 

546 final_noise_warps = { 

547 n_noise: self._prepareEmptyExposure(sky_info) 

548 for n_noise in range(self.config.numberOfNoiseRealizations) 

549 } 

550 

551 # We need a few bookkeeping variables only for the science coadd. 

552 totalGoodPixels = 0 

553 inputRecorder = self.inputRecorder.makeCoaddTempExpRecorder( 

554 visit_id, 

555 len(inputs), 

556 ) 

557 

558 for index, detector_inputs in enumerate(inputs.values()): 

559 self.log.debug( 

560 "Warping exposure %d/%d for id=%s", 

561 index + 1, 

562 len(inputs), 

563 detector_inputs.data_id, 

564 ) 

565 

566 input_exposure = detector_inputs.exposure 

567 # Generate noise image(s) in-situ. 

568 seed = self.get_seed_from_data_id(detector_inputs.data_id) 

569 # Limit to last 32 bits to avoid overflow in numpy. 

570 np_seed = (seed + self.config.seedOffset) & 0xFFFFFFFF 

571 self.log.debug("Setting numpy random seed to %d for noise realization", np_seed) 

572 rng = np.random.RandomState(np_seed) 

573 

574 # Generate noise images in-situ. 

575 noise_calexps = self.make_noise_exposures(input_exposure, rng) 

576 

577 warpedExposure = self.process( 

578 detector_inputs, 

579 target_wcs, 

580 self.warper, 

581 visit_summary, 

582 destBBox=target_bbox, 

583 ) 

584 

585 if warpedExposure is None: 

586 self.log.debug( 

587 "Skipping exposure %s because it could not be warped.", detector_inputs.data_id 

588 ) 

589 continue 

590 

591 # Accumulate the partial warps in an online fashion. 

592 nGood = copyGoodPixels( 

593 final_warp.maskedImage, 

594 warpedExposure.maskedImage, 

595 final_warp.mask.getPlaneBitMask(["NO_DATA"]), 

596 ) 

597 

598 ccdId = self.config.idGenerator.apply(detector_inputs.data_id).catalog_id 

599 inputRecorder.addCalExp(input_exposure, ccdId, nGood) 

600 totalGoodPixels += nGood 

601 

602 if self.config.doWarpMaskedFraction: 

603 # Obtain the masked fraction exposure and warp it. 

604 if self.config.doPreWarpInterpolation: 

605 badMaskPlanes = self.preWarpInterpolation.config.badMaskPlanes 

606 else: 

607 badMaskPlanes = [] 

608 masked_fraction_exp = self._get_bad_mask(input_exposure, badMaskPlanes) 

609 

610 masked_fraction_warp = self.maskedFractionWarper.warpExposure( 

611 target_wcs, masked_fraction_exp, destBBox=target_bbox 

612 ) 

613 

614 copyGoodPixels( 

615 final_masked_fraction_warp.maskedImage, 

616 masked_fraction_warp.maskedImage, 

617 final_masked_fraction_warp.mask.getPlaneBitMask(["NO_DATA"]), 

618 ) 

619 

620 # Process and accumulate noise images. 

621 for n_noise in range(self.config.numberOfNoiseRealizations): 

622 noise_calexp = noise_calexps[n_noise] 

623 noise_pseudo_inputs = dataclasses.replace( 

624 detector_inputs, 

625 exposure_or_handle=noise_calexp, 

626 background_revert=BackgroundList(), 

627 background_apply=BackgroundList(), 

628 ) 

629 warpedNoise = self.process( 

630 noise_pseudo_inputs, 

631 target_wcs, 

632 self.warper, 

633 visit_summary, 

634 destBBox=target_bbox, 

635 ) 

636 

637 copyGoodPixels( 

638 final_noise_warps[n_noise].maskedImage, 

639 warpedNoise.maskedImage, 

640 final_noise_warps[n_noise].mask.getPlaneBitMask(["NO_DATA"]), 

641 ) 

642 

643 # If there are no good pixels, return a Struct filled with None. 

644 if totalGoodPixels == 0: 

645 results = Struct( 

646 warp=None, 

647 masked_fraction_warp=None, 

648 ) 

649 for noise_index in range(self.config.numberOfNoiseRealizations): 

650 setattr(results, f"noise_warp{noise_index}", None) 

651 

652 return results 

653 

654 # Finish the inputRecorder and add the coaddPsf to the final warp. 

655 inputRecorder.finish(final_warp, totalGoodPixels) 

656 

657 coaddPsf = CoaddPsf( 

658 inputRecorder.coaddInputs.ccds, 

659 sky_info.wcs, 

660 self.config.coaddPsf.makeControl(), 

661 ) 

662 

663 final_warp.setPsf(coaddPsf) 

664 for _, warp_detector_input in inputs.items(): 

665 final_warp.setFilter(warp_detector_input.exposure.getFilter()) 

666 final_warp.getInfo().setVisitInfo(warp_detector_input.exposure.getInfo().getVisitInfo()) 

667 break # Just need the filter and visit info from any one of the inputs. 

668 

669 results = Struct( 

670 warp=final_warp, 

671 ) 

672 

673 if self.config.doWarpMaskedFraction: 

674 results.masked_fraction_warp = final_masked_fraction_warp.image 

675 

676 for noise_index, noise_exposure in final_noise_warps.items(): 

677 setattr(results, f"noise_warp{noise_index}", noise_exposure.maskedImage) 

678 

679 return results 

680 

681 def process( 

682 self, 

683 detector_inputs: WarpDetectorInputs, 

684 target_wcs, 

685 warper, 

686 visit_summary=None, 

687 maxBBox=None, 

688 destBBox=None, 

689 ) -> ExposureF | None: 

690 """Process an exposure. 

691 

692 There are three processing steps that are applied to the input: 

693 

694 1. Interpolate over bad pixels before warping. 

695 2. Apply all calibrations from visit_summary to the exposure. 

696 3. Warp the exposure to the target coordinate system. 

697 

698 Parameters 

699 ---------- 

700 detector_inputs : `WarpDetectorInputs` 

701 The input exposure to be processed, along with any other 

702 per-detector modifications. 

703 target_wcs : `~lsst.afw.geom.SkyWcs` 

704 The WCS of the target patch. 

705 warper : `~lsst.afw.math.Warper` 

706 The warper to use for warping the input exposure. 

707 visit_summary : `~lsst.afw.table.ExposureCatalog` | None 

708 Table of visit summary information. If not None, the visit_summary 

709 information will be used to update the calibration of the input 

710 exposures. Otherwise, the input exposures will be used as-is. 

711 maxBBox : `~lsst.geom.Box2I` | None 

712 Maximum bounding box of the warped exposure. If None, this is 

713 determined automatically. 

714 destBBox : `~lsst.geom.Box2I` | None 

715 Exact bounding box of the warped exposure. If None, this is 

716 determined automatically. 

717 

718 Returns 

719 ------- 

720 warped_exposure : `~lsst.afw.image.Exposure` | None 

721 The processed and warped exposure, if all the calibrations could 

722 be applied successfully. Otherwise, None. 

723 """ 

724 

725 input_exposure = detector_inputs.exposure 

726 if self.config.doPreWarpInterpolation: 

727 self.preWarpInterpolation.run(input_exposure.maskedImage) 

728 

729 success = self._apply_all_calibrations(detector_inputs, visit_summary=visit_summary) 

730 

731 if not success: 

732 return None 

733 

734 with self.timer("warp"): 

735 warped_exposure = warper.warpExposure( 

736 target_wcs, 

737 input_exposure, 

738 maxBBox=maxBBox, 

739 destBBox=destBBox, 

740 ) 

741 

742 # Potentially a post-warp interpolation here? Relies on DM-38630. 

743 

744 return warped_exposure 

745 

746 def _apply_all_calibrations( 

747 self, 

748 detector_inputs: WarpDetectorInputs, 

749 *, 

750 visit_summary: ExposureCatalog | None = None, 

751 ) -> bool: 

752 """Apply all of the calibrations from visit_summary to the exposure. 

753 

754 Specifically, this method updates the following (if available) to the 

755 input exposure in place from ``visit_summary``: 

756 

757 - Aperture correction map 

758 - Photometric calibration 

759 - PSF 

760 - WCS 

761 

762 It also reverts and applies backgrounds in ``detector_inputs``. 

763 

764 Parameters 

765 ---------- 

766 detector_inputs : `WarpDetectorInputs` 

767 The input exposure to be processed, along with any other 

768 per-detector modifications. 

769 visit_summary : `~lsst.afw.table.ExposureCatalog` | None 

770 Table of visit summary information. If not None, the visit summary 

771 information will be used to update the calibration of the input 

772 exposures. Otherwise, the input exposures will be used as-is. 

773 

774 Returns 

775 ------- 

776 success : `bool` 

777 True if all calibrations were successfully applied, 

778 False otherwise. 

779 

780 Raises 

781 ------ 

782 RuntimeError 

783 Raised if ``visit_summary`` is provided but does not contain a 

784 record corresponding to ``detector_inputs``, or if one of the 

785 background fields in ``detector_inputs`` is inconsistent with the 

786 task configuration. 

787 """ 

788 if self.config.doRevertOldBackground: 

789 detector_inputs.revert_background() 

790 elif detector_inputs.background_revert: 

791 # This could get trigged only if runQuantum is skipped and run is 

792 # called directly. 

793 raise RuntimeError( 

794 f"doRevertOldBackground is False, but {detector_inputs.data_id} has a background_revert." 

795 ) 

796 

797 input_exposure = detector_inputs.exposure 

798 if visit_summary is not None: 

799 detector = input_exposure.info.getDetector().getId() 

800 row = visit_summary.find(detector) 

801 

802 if row is None: 

803 self.log.info( 

804 "Unexpectedly incomplete visit_summary: detector = %s is missing. Skipping it.", 

805 detector, 

806 ) 

807 return False 

808 

809 if self.config.useVisitSummaryPhotoCalib: 

810 photo_calib = row.getPhotoCalib() 

811 if photo_calib is None: 

812 self.log.info( 

813 "No photometric calibration found in visit summary for detector = %s. Skipping it.", 

814 detector, 

815 ) 

816 return False 

817 elif not np.isfinite(photo_calib.getInstFluxAtZeroMagnitude()): 

818 self.log.info( 

819 f"Invalid PhotoCalib found in visit summary for detector {detector}. Skipping it.", 

820 ) 

821 return False 

822 else: 

823 input_exposure.setPhotoCalib(photo_calib) 

824 

825 if self.config.useVisitSummaryWcs: 

826 if wcs := row.getWcs(): 

827 input_exposure.setWcs(wcs) 

828 else: 

829 self.log.info("No WCS found in visit summary for detector = %s. Skipping it.", detector) 

830 return False 

831 

832 if self.config.useVisitSummaryPsf: 

833 if psf := row.getPsf(): 

834 input_exposure.setPsf(psf) 

835 else: 

836 self.log.info("No PSF found in visit summary for detector = %s. Skipping it.", detector) 

837 return False 

838 

839 if apcorr_map := row.getApCorrMap(): 

840 input_exposure.setApCorrMap(apcorr_map) 

841 else: 

842 self.log.info( 

843 "No aperture correction map found in visit summary for detector = %s. Skipping it", 

844 detector, 

845 ) 

846 return False 

847 

848 if self.config.doApplyNewBackground: 

849 detector_inputs.apply_background() 

850 elif detector_inputs.background_apply: 

851 # This could get trigged only if runQuantum is skipped and run is 

852 # called directly. 

853 raise RuntimeError( 

854 f"doApplyNewBackground is False, but {detector_inputs.data_id} has a background_apply." 

855 ) 

856 

857 # Calibrate the (masked) image. 

858 # This should likely happen even if visit_summary is None. 

859 photo_calib = input_exposure.photoCalib 

860 input_exposure.maskedImage = photo_calib.calibrateImage(input_exposure.maskedImage) 

861 # This new PhotoCalib shouldn't need to be used, but setting it here 

862 # to reflect the fact that the image now has calibrated pixels might 

863 # help avoid future bugs. 

864 input_exposure.setPhotoCalib(PhotoCalib(1.0)) 

865 

866 return True 

867 

868 # This method is copied from makeWarp.py 

869 @classmethod 

870 def _prepareEmptyExposure(cls, sky_info): 

871 """Produce an empty exposure for a given patch. 

872 

873 Parameters 

874 ---------- 

875 sky_info : `lsst.pipe.base.Struct` 

876 Struct from `~lsst.pipe.base.coaddBase.makeSkyInfo` with 

877 geometric information about the patch. 

878 

879 Returns 

880 ------- 

881 exp : `lsst.afw.image.exposure.ExposureF` 

882 An empty exposure for a given patch. 

883 """ 

884 exp = ExposureF(sky_info.bbox, sky_info.wcs) 

885 exp.getMaskedImage().set(np.nan, Mask.getPlaneBitMask("NO_DATA"), np.inf) 

886 # Set the PhotoCalib to 1 to mean that pixels are nJy, since we will 

887 # calibrate them before we warp them. 

888 exp.setPhotoCalib(PhotoCalib(1.0)) 

889 exp.metadata["BUNIT"] = "nJy" 

890 return exp 

891 

892 @staticmethod 

893 def compute_median_variance(mi: MaskedImage) -> float: 

894 """Compute the median variance across the good pixels of a MaskedImage. 

895 

896 Parameters 

897 ---------- 

898 mi : `~lsst.afw.image.MaskedImage` 

899 The input image on which to compute the median variance. 

900 

901 Returns 

902 ------- 

903 median_variance : `float` 

904 Median variance of the input calexp. 

905 """ 

906 # Shouldn't this exclude pixels that are masked, to be safe? 

907 # This is implemented as it was in descwl_coadd. 

908 return np.median(mi.variance.array[np.isfinite(mi.variance.array) & np.isfinite(mi.image.array)]) 

909 

910 def get_seed_from_data_id(self, data_id) -> int: 

911 """Get a seed value given a data_id. 

912 

913 This method generates a unique, reproducible pseudo-random number for 

914 a data id. This is not affected by ordering of the input, or what 

915 set of visits, ccds etc. are given. 

916 

917 This is implemented as a public method, so that simulations that 

918 don't necessary deal with the middleware can mock up a ``data_id`` 

919 instance, or override this method with a different one to obtain a 

920 seed value consistent with the pipeline task. 

921 

922 Parameters 

923 ---------- 

924 data_id : `~lsst.daf.butler.DataCoordinate` 

925 Data identifier dictionary. 

926 

927 Returns 

928 ------- 

929 seed : `int` 

930 A unique seed for this data_id to seed a random number generator.