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

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 astropy.units 

25from astropy.stats import gaussian_sigma_to_fwhm 

26import numpy as np 

27 

28import lsst.afw.detection as afwDetection 

29import lsst.afw.image as afwImage 

30import lsst.geom as geom 

31import lsst.afw.geom as afwGeom 

32from lsst.afw.image import VisitInfo 

33import lsst.afw.table as afwTable 

34from lsst.afw.math._warper import computeWarpedBBox 

35import lsst.afw.math as afwMath 

36import lsst.pex.config as pexConfig 

37import lsst.pex.exceptions 

38import lsst.pipe.base as pipeBase 

39 

40from lsst.images import ( 

41 DifferenceImage, 

42 DifferenceImageTemplateInfo, 

43 get_legacy_optional_mask_planes, 

44 get_legacy_template_mask_planes, 

45) 

46from lsst.skymap import BaseSkyMap 

47from lsst.ip.diffim.dcrModel import DcrModel 

48from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm, record_from_obs_info 

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

50from lsst.utils.timer import timeMethod 

51 

52__all__ = [ 

53 "GetTemplateTask", 

54 "GetTemplateConfig", 

55 "GetDcrTemplateTask", 

56 "GetDcrTemplateConfig", 

57] 

58 

59_PSF_STAMP_SIGMAS = 5 

60"""Width, in N sigma, of the Gaussian approximation to the CoaddPsf.""" 

61 

62 

63class GetTemplateConnections( 

64 pipeBase.PipelineTaskConnections, 

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

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

67): 

68 bbox = pipeBase.connectionTypes.Input( 

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

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

71 storageClass="Box2I", 

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

73 deprecated="Replaced by the `detector` connection. Will be removed after v31.", 

74 ) 

75 detector = pipeBase.connectionTypes.Input( 

76 doc="Detector of the exposure that we will construct the template for." 

77 " Its bounding box sets the geometry of the output template.", 

78 name="{fakesType}calexp.detector", 

79 storageClass="Detector", 

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

81 ) 

82 obs_info = pipeBase.connectionTypes.Input( 

83 doc="Observation metadata of the exposure the template is built for." 

84 " Only read with output_image_type='future', which needs this to set" 

85 " the obs_info component of a lsst.images.VisitImage.", 

86 name="{fakesType}calexp.obs_info", 

87 storageClass="ObservationInfo", 

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

89 ) 

90 wcs = pipeBase.connectionTypes.Input( 

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

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

93 storageClass="Wcs", 

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

95 ) 

96 skyMap = pipeBase.connectionTypes.Input( 

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

98 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

99 dimensions=("skymap",), 

100 storageClass="SkyMap", 

101 ) 

102 coaddExposures = pipeBase.connectionTypes.Input( 

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

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

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

106 storageClass="ExposureF", 

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

108 multiple=True, 

109 deferLoad=True, 

110 deferGraphConstraint=True, 

111 ) 

112 

113 template = pipeBase.connectionTypes.Output( 

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

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

116 storageClass="ExposureF", 

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

118 ) 

119 

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

121 super().__init__(config=config) 

122 # Kept only so that older config files still load. 

123 del self.bbox 

124 if config.requireCoaddAtGraphBuild: 124 ↛ 125line 124 didn't jump to line 125 because the condition on line 124 was never true

125 self.coaddExposures = dataclasses.replace( 

126 self.coaddExposures, 

127 deferGraphConstraint=False, 

128 ) 

129 if config.output_image_type == "future": 

130 self.template = dataclasses.replace(self.template, storageClass="DifferenceImage") 

131 else: 

132 del self.obs_info 

133 

134 

135class GetTemplateConfig( 

136 pipeBase.PipelineTaskConfig, pipelineConnections=GetTemplateConnections 

137): 

138 templateBorderSize = pexConfig.Field( 

139 dtype=int, 

140 default=20, 

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

142 ) 

143 warp = pexConfig.ConfigField( 

144 dtype=afwMath.Warper.ConfigClass, 

145 doc="warper configuration", 

146 ) 

147 doCorrectVariancePlateScale = pexConfig.Field( 

148 dtype=bool, 

149 default=False, 

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

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

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

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

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

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

156 ) 

157 coaddPsf = pexConfig.ConfigField( 

158 doc="Configuration for CoaddPsf", 

159 dtype=CoaddPsfConfig, 

160 ) 

161 fwhmExposureGrid = pexConfig.Field( 

162 doc="Grid size to compute the average PSF FWHM in an exposure (pixels)", 

163 dtype=int, 

164 default=10, 

165 ) 

166 fwhmExposureBuffer = pexConfig.Field( 

167 doc="Fractional buffer margin to be left out of all sides of the image during construction" 

168 " of grid to compute average PSF FWHM in an exposure", 

169 dtype=float, 

170 default=0.05, 

171 ) 

172 varianceBackground = pexConfig.ConfigurableField( 

173 target=SubtractBackgroundTask, 

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

175 ) 

176 highVarianceThreshold = pexConfig.RangeField( 

177 dtype=float, 

178 default=4, 

179 min=1, 

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

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

182 ) 

183 highVarianceMaskFraction = pexConfig.Field( 

184 dtype=float, 

185 default=0.1, 

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

187 " HIGH_VARIANCE mask plane.", 

188 ) 

189 doScaleVariance = pexConfig.Field( 

190 dtype=bool, 

191 default=True, 

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

193 ) 

194 scaleVariance = pexConfig.ConfigurableField( 

195 target=ScaleVarianceTask, 

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

197 ) 

198 requireCoaddAtGraphBuild = pexConfig.Field( 

199 dtype=bool, 

200 default=False, 

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

202 " initial butler query during QuantumGraph generation.", 

203 ) 

204 output_image_type = pexConfig.ChoiceField[str]( 

205 "Which image type to use for the output template.", 

206 allowed={ 

207 "legacy": "Write as a lsst.afw.image.ExposureF.", 

208 "future": "Write as a lsst.images.DifferenceImage.", 

209 }, 

210 optional=False, 

211 default="legacy", 

212 ) 

213 raiseOnUndefinedMaskMap = pexConfig.Field( 

214 dtype=bool, 

215 default=False, 

216 doc="Raise if the template has pixels set in a mask plane that" 

217 " `lsst.images` does not map? If False, log a warning and drop that" 

218 " plane from the output instead." 

219 " Only used with output_image_type='future'.", 

220 ) 

221 

222 def setDefaults(self): 

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

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

225 self.warp.cacheSize = 100000 

226 self.coaddPsf.cacheSize = self.warp.cacheSize 

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

228 # interpolation length for WCS evaluations. 

229 self.warp.interpLength = 100 

230 self.warp.warpingKernelName = "lanczos3" 

231 self.coaddPsf.warpingKernelName = self.warp.warpingKernelName 

232 

233 # Background subtraction of the variance plane 

234 self.varianceBackground.algorithm = "LINEAR" 

235 self.varianceBackground.binSize = 32 

236 self.varianceBackground.useApprox = False 

237 self.varianceBackground.statisticsProperty = "MEDIAN" 

238 self.varianceBackground.doFilterSuperPixels = True 

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

240 "EDGE", 

241 "DETECTED", 

242 "DETECTED_NEGATIVE", 

243 "NO_DATA", 

244 ] 

245 

246 

247class GetTemplateTask(pipeBase.PipelineTask): 

248 ConfigClass = GetTemplateConfig 

249 _DefaultName = "getTemplate" 

250 

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

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

253 if self.config.doScaleVariance: 

254 self.makeSubtask("scaleVariance") 

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

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

257 self.schema.addField( 

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

259 ) 

260 self.schema.addField( 

261 "patch", 

262 type=np.int32, 

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

264 ) 

265 self.schema.addField( 

266 "weight", 

267 type=float, 

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

269 ) 

270 self.makeSubtask("varianceBackground") 

271 

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

273 inputs = butlerQC.get(inputRefs) 

274 detector = inputs.pop("detector") 

275 bbox = detector.getBBox() 

276 wcs = inputs.pop("wcs") 

277 obsInfo = inputs.pop("obs_info", None) 

278 coaddExposures = inputs.pop("coaddExposures") 

279 skymap = inputs.pop("skyMap") 

280 

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

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

283 

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

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

286 outputs = self.run( 

287 coaddExposureHandles=results.coaddExposures, 

288 bbox=bbox, 

289 wcs=wcs, 

290 dataIds=results.dataIds, 

291 physical_filter=physical_filter, 

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

293 ) 

294 if self.config.output_image_type == "future": 

295 dataId = butlerQC.quantum.dataId 

296 exposureRecord = record_from_obs_info(obsInfo, dataId["instrument"], dataId["visit"], 

297 dataId.universe) 

298 outputs.template = self.convert_outputs_to_future(outputs.template, 

299 inputRefs.coaddExposures, 

300 detector=detector, 

301 exposureRecord=exposureRecord) 

302 butlerQC.put(outputs, outputRefs) 

303 

304 def convert_outputs_to_future(self, template, coadd_refs, detector, exposureRecord): 

305 """Convert a template image to DifferenceImage format. 

306 

307 This replaces ``template`` with an `lsst.images.DifferenceImage` 

308 instance and records the coadds contributed to it. The Gaussian PSF 

309 set by `run` is converted to an 

310 `lsst.images.psfs.GaussianPointSpreadFunction`. 

311 

312 Parameters 

313 ---------- 

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

315 Template to convert. Its detector is replaced with ``detector`` 

316 and its mask is modified in place. 

317 coadd_refs : `list` [`lsst.daf.butler.DatasetRef`] 

318 References to the coadds that may have gone into the template. 

319 May be a superset of the coadds that were used; supplies the 

320 dataset id and RUN collection recorded for each coadd that was. 

321 detector : `lsst.afw.cameraGeom.Detector` 

322 Detector of the science image that the template was built for. 

323 exposureRecord : `lsst.daf.butler.DimensionRecord` 

324 The ``exposure`` record of the observation the template was built 

325 for, which supplies the observation metadata recorded on the 

326 output. 

327 

328 Returns 

329 ------- 

330 convertedTemplate : `lsst.images.DifferenceImage` 

331 The converted image. 

332 

333 Raises 

334 ------ 

335 RuntimeError 

336 Raised if the template has pixels set in a mask plane that 

337 `lsst.images.get_legacy_template_mask_planes` does not map and 

338 ``config.raiseOnUndefinedMaskMap`` is True, or if ``coadd_refs`` 

339 spans more than one skymap. 

340 """ 

341 # The template PSF is a Gaussian, so rebuild the CoaddPsf from the 

342 # template inputs to record which coadds contributed. 

343 coaddPsf = self._makePsf(template, template.getInfo().getCoaddInputs().ccds, template.getWcs()) 

344 template.setDetector(detector) 

345 skymaps = {ref.dataId["skymap"] for ref in coadd_refs} 

346 if len(skymaps) != 1: 346 ↛ 347line 346 didn't jump to line 347 because the condition on line 346 was never true

347 raise RuntimeError("Expected the template inputs to come from one skymap;" 

348 f" got {sorted(skymaps)}.") 

349 butler_info = {(ref.dataId["tract"], ref.dataId["patch"]): (ref.id, ref.run) 

350 for ref in coadd_refs} 

351 planeMap = get_legacy_template_mask_planes() 

352 if not self.config.raiseOnUndefinedMaskMap: 

353 self._clearUnmappedMaskPlanes(template.mask, planeMap) 

354 convertedTemplate = DifferenceImage.from_legacy( 

355 template, 

356 exposure_record=exposureRecord, 

357 unit=astropy.units.nJy, 

358 plane_map=planeMap, 

359 ) 

360 del template 

361 convertedTemplate.templates = DifferenceImageTemplateInfo.from_legacy_psf( 

362 convertedTemplate.sky_projection.pixel_frame, 

363 coaddPsf, 

364 skymaps.pop(), 

365 butler_info, 

366 log=self.log, 

367 ) 

368 return convertedTemplate 

369 

370 def _clearUnmappedMaskPlanes(self, mask, planeMap): 

371 """Clear the pixels of every mask plane that conversion would reject. 

372 

373 Parameters 

374 ---------- 

375 mask : `lsst.afw.image.Mask` 

376 Template mask to modify in place. 

377 planeMap : `dict` [`str`, `lsst.images.MaskPlane`] 

378 Plane map the template will be converted with. The optional 

379 planes of `lsst.images.get_legacy_optional_mask_planes` are kept 

380 as well, since conversion adds them when they have pixels set. 

381 """ 

382 keep = planeMap.keys() | get_legacy_optional_mask_planes().keys() 

383 for name in mask.getMaskPlaneDict(): 

384 if name in keep: 

385 continue 

386 bitmask = mask.getPlaneBitMask(name) 

387 nSet = np.count_nonzero(mask.array & bitmask) 

388 if nSet: 

389 self.log.warning("Dropping mask plane %s from the template: %d pixels are set, but" 

390 " the template plane map does not include it.", name, nSet) 

391 mask.array &= ~bitmask 

392 

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

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

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

396 structure of their dataIds. 

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

398 

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

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

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

402 all possibly-matching coadd exposures. 

403 

404 Parameters 

405 ---------- 

406 coaddExposureHandles : `iterable` \ 

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

408 `lsst.afw.image.Exposure`] 

409 Dataset handles to exposures that might overlap the desired 

410 region. 

411 bbox : `lsst.geom.Box2I` 

412 Template bounding box of the pixel geometry onto which the 

413 coaddExposures will be resampled. 

414 skymap : `lsst.skymap.SkyMap` 

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

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

417 Template WCS onto which the coadds will be resampled. 

418 

419 Returns 

420 ------- 

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

422 A struct with attributes: 

423 

424 ``coaddExposures`` 

425 Dict of coadd exposures that overlap the projected bbox, 

426 indexed on tract id 

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

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

429 ``dataIds`` 

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

431 projected bbox, indexed on tract id 

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

433 

434 Raises 

435 ------ 

436 NoWorkFound 

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

438 WCS is None. 

439 """ 

440 if wcs is None: 

441 raise pipeBase.NoWorkFound( 

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

443 ) 

444 

445 # Exposure's validPolygon would be more accurate 

446 detectorPolygon = geom.Box2D(bbox) 

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

448 overlappingArea = 0 

449 coaddExposures = collections.defaultdict(list) 

450 dataIds = collections.defaultdict(list) 

451 

452 for coaddRef in coaddExposureHandles: 

453 dataId = coaddRef.dataId 

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

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

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

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

458 # detector coordinates because the skymap's inverse mapping 

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

460 # arbitrary sky coordinates. 

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

462 if patchPolygon.intersection(detectorInPatchCoordinates): 

463 overlappingArea += patchPolygon.intersectionSingle( 

464 detectorInPatchCoordinates 

465 ).calculateArea() 

466 self.log.info( 

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

468 dataId["tract"], 

469 dataId["patch"], 

470 ) 

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

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

473 

474 if not overlappingArea: 

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

476 

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

478 

479 @timeMethod 

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

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

482 subtract from a science image. 

483 

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

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

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

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

488 The PSF of the template is a Gaussian with the same width as the 

489 CoaddPsf evaluated at the center of the detector. 

490 

491 Parameters 

492 ---------- 

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

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

495 `lsst.afw.image.Exposure`]] 

496 Coadds to be mosaicked, indexed on tract id. 

497 bbox : `lsst.geom.Box2I` 

498 Input Bounding box of the detector geometry onto which to 

499 resample the ``coaddExposureHandles``. The output template bbox 

500 will be grown by ``templateBorderSize``. 

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

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

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

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

505 tract id. 

506 physical_filter : `str` 

507 Physical filter of the science image. 

508 visit : `int`, optional 

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

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

511 science image for the visit. 

512 

513 Returns 

514 ------- 

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

516 A struct with attributes: 

517 

518 ``template`` 

519 A template coadd exposure assembled out of patches, with 

520 pixels in nJy. Its PSF is a Gaussian approximation to the 

521 CoaddPsf of its inputs (`lsst.afw.image.ExposureF`). 

522 

523 Raises 

524 ------ 

525 NoWorkFound 

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

527 """ 

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

529 

530 bbox = geom.Box2I(bbox) 

531 bbox.grow(self.config.templateBorderSize) 

532 

533 warped = {} 

534 catalogs = [] 

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

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

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

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

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

540 # different instruments are not supported. 

541 plateScaleFactor = None 

542 if self.config.doCorrectVariancePlateScale and coaddExposureHandles: 

543 plateScaleFactor = self._plateScaleFactor( 

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

545 ) 

546 

547 for tract in coaddExposureHandles: 

548 maskedImages, catalog, totalBox = self._makeExposureCatalog( 

549 coaddExposureHandles[tract], dataIds[tract] 

550 ) 

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

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

553 # Combine images from individual patches together. 

554 unwarped, count, included = self._merge( 

555 maskedImages, warpedBox, catalog[0].wcs 

556 ) 

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

558 del maskedImages 

559 if count == 0: 

560 self.log.info( 

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

562 tract, 

563 ) 

564 continue 

565 warpedBox.clip(totalBox) 

566 potentialInput = self.warper.warpExposure( 

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

568 ) 

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

570 del unwarped 

571 if np.all( 

572 potentialInput.mask.array 

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

574 ): 

575 self.log.info( 

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

577 tract, 

578 ) 

579 continue 

580 

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

582 tempCatalog = afwTable.ExposureCatalog(self.schema) 

583 tempCatalog.reserve(len(included)) 

584 for i in included: 

585 tempCatalog.append(catalog[i]) 

586 catalogs.append(tempCatalog) 

587 warped[tract] = potentialInput.maskedImage 

588 

589 if len(warped) == 0: 

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

591 

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

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

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

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

596 

597 if self.config.doScaleVariance: 

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

599 # needed. Note that the science variance is scaled 

600 # independently in ``AlardLuptonSubtractTask``. 

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

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

603 self.metadata["scaleTemplateVarianceFactor"] = varianceFactor 

604 

605 self._correctVariance(template, plateScaleFactor) 

606 

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

608 catalog = afwTable.ExposureCatalog(self.schema) 

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

610 for c in catalogs: 

611 catalog.extend(c) 

612 

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

614 self.checkHighVariance(template) 

615 if visit is not None: 

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

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

618 identity = afwImage.PhotoCalib(1.0) 

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

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

621 photoCalib = identity 

622 template.setPhotoCalib(photoCalib) 

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

624 template.setPsf(self._makeGaussianPsf(template)) 

625 

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

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

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

629 template.getInfo().setCoaddInputs(coaddInputs) 

630 return pipeBase.Struct(template=template) 

631 

632 def _correctVariance(self, template, plateScaleFactor): 

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

634 

635 Compute a multiplicative correction factor for the template variance 

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

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

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

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

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

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

642 image pixel noise level would suggest. 

643 

644 Parameters 

645 ---------- 

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

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

648 plateScaleFactor : `float` or `None` 

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

650 

651 Raises 

652 ------ 

653 RuntimeError 

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

655 could not be reconstructed from the coadd inputs. 

656 """ 

657 scale = 1.0 

658 if self.config.doCorrectVariancePlateScale: 

659 if plateScaleFactor is None: 

660 raise RuntimeError( 

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

662 ) 

663 scale *= plateScaleFactor 

664 self.metadata["variancePlateScaleFactor"] = plateScaleFactor 

665 self.log.info( 

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

667 " coadd inputs.", 

668 plateScaleFactor, 

669 ) 

670 

671 if scale == 1.0: 

672 return 

673 template.variance.array *= scale 

674 self.metadata["templateVarianceCorrectionFactor"] = scale 

675 self.log.info( 

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

677 scale, 

678 ) 

679 

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

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

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

683 

684 Parameters 

685 ---------- 

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

687 `lsst.afw.image.Exposure` 

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

689 component is read; the pixels are left alone. 

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

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

692 bbox : `lsst.geom.Box2I` 

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

694 evaluate the science image pixel scale. 

695 

696 Returns 

697 ------- 

698 factor : `float` or `None` 

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

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

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

702 """ 

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

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

705 return None 

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

707 

708 for record in coaddInputs.ccds: 

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

710 # that gives a usable pixel area. 

711 originalWcs = record.getWcs() 

712 if originalWcs is None: 

713 continue 

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

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

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

717 factor = (pviScale/scienceScale)**2 

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

719 return factor 

720 return None 

721 

722 def checkHighVariance(self, template): 

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

724 

725 Parameters 

726 ---------- 

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

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

729 """ 

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

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

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

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

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

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

736 " pixels were unmasked for background estimation, but" 

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

738 else: 

739 varianceExposure = template.clone() 

740 varianceExposure.image.array = varianceExposure.variance.array 

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

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

743 highVariancePix = varianceBackground > threshold 

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

745 

746 @staticmethod 

747 def _checkInputs(dataIds, coaddExposures): 

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

749 the exposures all have the same photometric calibration. 

750 

751 Parameters 

752 ---------- 

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

754 Record of the tract and patch of each coaddExposure. 

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

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

757 `lsst.afw.image.Exposure` or 

758 `lsst.afw.image.Exposure`]] 

759 Coadds to be mosaicked. 

760 

761 Returns 

762 ------- 

763 band : `str` 

764 Filter band of all the input exposures. 

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

766 Photometric calibration of all of the input exposures. 

767 

768 Raises 

769 ------ 

770 RuntimeError 

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

772 all the same. 

773 """ 

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

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

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

777 band = bands.pop() 

778 photoCalibs = [ 

779 exposure.get(component="photoCalib") 

780 for exposures in coaddExposures.values() 

781 for exposure in exposures 

782 ] 

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

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

785 raise RuntimeError(msg) 

786 photoCalib = photoCalibs[0] 

787 return band, photoCalib 

788 

789 def _makeExposureCatalog(self, exposureRefs, dataIds): 

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

791 

792 Parameters 

793 ---------- 

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

795 `lsst.afw.image.Exposure`] 

796 Exposures to include in the catalog. 

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

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

799 "patch" entries. 

800 

801 Returns 

802 ------- 

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

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

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

806 Catalog of metadata for each exposure 

807 totalBox : `lsst.geom.Box2I` 

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

809 """ 

810 catalog = afwTable.ExposureCatalog(self.schema) 

811 catalog.reserve(len(exposureRefs)) 

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

813 images = {} 

814 totalBox = geom.Box2I() 

815 

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

817 images[dataId] = coadd.maskedImage 

818 bbox = coadd.getBBox() 

819 totalBox = totalBox.expandedTo(bbox) 

820 record = catalog.addNew() 

821 record.setPsf(coadd.psf) 

822 record.setWcs(coadd.wcs) 

823 record.setPhotoCalib(coadd.photoCalib) 

824 record.setBBox(bbox) 

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

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

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

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

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

830 record.set("weight", 1) 

831 

832 return images, catalog, totalBox 

833 

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

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

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

837 exposures. 

838 

839 Parameters 

840 ---------- 

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

842 `lsst.afw.image.Exposure`] 

843 Images to be merged into one larger bounding box. 

844 bbox : `lsst.geom.Box2I` 

845 Bounding box defining the image to merge into. 

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

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

848 

849 Returns 

850 ------- 

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

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

853 positions. 

854 count : `int` 

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

856 in the merged image. 

857 included : `list` [`int`] 

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

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

860 """ 

861 merged = afwImage.ExposureF(bbox, wcs) 

862 weights = afwImage.ImageF(bbox) 

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

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

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

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

867 clippedBox.clip(bbox) 

868 if clippedBox.area == 0: 

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

870 continue # nothing in this image overlaps the output 

871 maskedImage = maskedImage.subset(clippedBox) 

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

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

874 np.isfinite(maskedImage.variance.array) 

875 ) 

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

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

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

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

880 maskedImage.image.array[bad] = 0.0 

881 maskedImage.variance.array[bad] = 0.0 

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

883 maskedImage.mask.array[bad] = 0 

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

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

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

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

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

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