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

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 

24from astropy.stats import gaussian_sigma_to_fwhm 

25import numpy as np 

26 

27import lsst.afw.detection as afwDetection 

28import lsst.afw.image as afwImage 

29import lsst.geom as geom 

30import lsst.afw.geom as afwGeom 

31from lsst.afw.image import VisitInfo 

32import lsst.afw.table as afwTable 

33from lsst.afw.math._warper import computeWarpedBBox 

34import lsst.afw.math as afwMath 

35import lsst.pex.config as pexConfig 

36import lsst.pex.exceptions 

37import lsst.pipe.base as pipeBase 

38 

39from lsst.skymap import BaseSkyMap 

40from lsst.ip.diffim.dcrModel import DcrModel 

41from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm 

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

43from lsst.utils.timer import timeMethod 

44 

45__all__ = [ 

46 "GetTemplateTask", 

47 "GetTemplateConfig", 

48 "GetDcrTemplateTask", 

49 "GetDcrTemplateConfig", 

50] 

51 

52_PSF_STAMP_SIGMAS = 5 

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

54 

55 

56class GetTemplateConnections( 

57 pipeBase.PipelineTaskConnections, 

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

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

60): 

61 bbox = pipeBase.connectionTypes.Input( 

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

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

64 storageClass="Box2I", 

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

66 ) 

67 wcs = pipeBase.connectionTypes.Input( 

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

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

70 storageClass="Wcs", 

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

72 ) 

73 skyMap = pipeBase.connectionTypes.Input( 

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

75 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

76 dimensions=("skymap",), 

77 storageClass="SkyMap", 

78 ) 

79 coaddExposures = pipeBase.connectionTypes.Input( 

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

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

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

83 storageClass="ExposureF", 

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

85 multiple=True, 

86 deferLoad=True, 

87 deferGraphConstraint=True, 

88 ) 

89 

90 template = pipeBase.connectionTypes.Output( 

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

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

93 storageClass="ExposureF", 

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

95 ) 

96 

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

98 super().__init__(config=config) 

99 if config.requireCoaddAtGraphBuild: 

100 self.coaddExposures = dataclasses.replace( 

101 self.coaddExposures, 

102 deferGraphConstraint=False, 

103 ) 

104 

105 

106class GetTemplateConfig( 

107 pipeBase.PipelineTaskConfig, pipelineConnections=GetTemplateConnections 

108): 

109 templateBorderSize = pexConfig.Field( 

110 dtype=int, 

111 default=20, 

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

113 ) 

114 warp = pexConfig.ConfigField( 

115 dtype=afwMath.Warper.ConfigClass, 

116 doc="warper configuration", 

117 ) 

118 doCorrectVariancePlateScale = pexConfig.Field( 

119 dtype=bool, 

120 default=False, 

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

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

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

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

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

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

127 ) 

128 coaddPsf = pexConfig.ConfigField( 

129 doc="Configuration for CoaddPsf", 

130 dtype=CoaddPsfConfig, 

131 ) 

132 fwhmExposureGrid = pexConfig.Field( 

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

134 dtype=int, 

135 default=10, 

136 ) 

137 fwhmExposureBuffer = pexConfig.Field( 

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

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

140 dtype=float, 

141 default=0.05, 

142 ) 

143 varianceBackground = pexConfig.ConfigurableField( 

144 target=SubtractBackgroundTask, 

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

146 ) 

147 highVarianceThreshold = pexConfig.RangeField( 

148 dtype=float, 

149 default=4, 

150 min=1, 

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

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

153 ) 

154 highVarianceMaskFraction = pexConfig.Field( 

155 dtype=float, 

156 default=0.1, 

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

158 " HIGH_VARIANCE mask plane.", 

159 ) 

160 doScaleVariance = pexConfig.Field( 

161 dtype=bool, 

162 default=True, 

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

164 ) 

165 scaleVariance = pexConfig.ConfigurableField( 

166 target=ScaleVarianceTask, 

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

168 ) 

169 requireCoaddAtGraphBuild = pexConfig.Field( 

170 dtype=bool, 

171 default=False, 

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

173 " initial butler query during QuantumGraph generation.", 

174 ) 

175 

176 def setDefaults(self): 

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

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

179 self.warp.cacheSize = 100000 

180 self.coaddPsf.cacheSize = self.warp.cacheSize 

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

182 # interpolation length for WCS evaluations. 

183 self.warp.interpLength = 100 

184 self.warp.warpingKernelName = "lanczos3" 

185 self.coaddPsf.warpingKernelName = self.warp.warpingKernelName 

186 

187 # Background subtraction of the variance plane 

188 self.varianceBackground.algorithm = "LINEAR" 

189 self.varianceBackground.binSize = 32 

190 self.varianceBackground.useApprox = False 

191 self.varianceBackground.statisticsProperty = "MEDIAN" 

192 self.varianceBackground.doFilterSuperPixels = True 

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

194 "EDGE", 

195 "DETECTED", 

196 "DETECTED_NEGATIVE", 

197 "NO_DATA", 

198 ] 

199 

200 

201class GetTemplateTask(pipeBase.PipelineTask): 

202 ConfigClass = GetTemplateConfig 

203 _DefaultName = "getTemplate" 

204 

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

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

207 if self.config.doScaleVariance: 

208 self.makeSubtask("scaleVariance") 

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

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

211 self.schema.addField( 

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

213 ) 

214 self.schema.addField( 

215 "patch", 

216 type=np.int32, 

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

218 ) 

219 self.schema.addField( 

220 "weight", 

221 type=float, 

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

223 ) 

224 self.makeSubtask("varianceBackground") 

225 

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

227 inputs = butlerQC.get(inputRefs) 

228 bbox = inputs.pop("bbox") 

229 wcs = inputs.pop("wcs") 

230 coaddExposures = inputs.pop("coaddExposures") 

231 skymap = inputs.pop("skyMap") 

232 

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

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

235 

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

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

238 outputs = self.run( 

239 coaddExposureHandles=results.coaddExposures, 

240 bbox=bbox, 

241 wcs=wcs, 

242 dataIds=results.dataIds, 

243 physical_filter=physical_filter, 

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

245 ) 

246 butlerQC.put(outputs, outputRefs) 

247 

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

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

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

251 structure of their dataIds. 

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

253 

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

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

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

257 all possibly-matching coadd exposures. 

258 

259 Parameters 

260 ---------- 

261 coaddExposureHandles : `iterable` \ 

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

263 `lsst.afw.image.Exposure`] 

264 Dataset handles to exposures that might overlap the desired 

265 region. 

266 bbox : `lsst.geom.Box2I` 

267 Template bounding box of the pixel geometry onto which the 

268 coaddExposures will be resampled. 

269 skymap : `lsst.skymap.SkyMap` 

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

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

272 Template WCS onto which the coadds will be resampled. 

273 

274 Returns 

275 ------- 

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

277 A struct with attributes: 

278 

279 ``coaddExposures`` 

280 Dict of coadd exposures that overlap the projected bbox, 

281 indexed on tract id 

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

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

284 ``dataIds`` 

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

286 projected bbox, indexed on tract id 

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

288 

289 Raises 

290 ------ 

291 NoWorkFound 

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

293 WCS is None. 

294 """ 

295 if wcs is None: 

296 raise pipeBase.NoWorkFound( 

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

298 ) 

299 

300 # Exposure's validPolygon would be more accurate 

301 detectorPolygon = geom.Box2D(bbox) 

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

303 overlappingArea = 0 

304 coaddExposures = collections.defaultdict(list) 

305 dataIds = collections.defaultdict(list) 

306 

307 for coaddRef in coaddExposureHandles: 

308 dataId = coaddRef.dataId 

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

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

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

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

313 # detector coordinates because the skymap's inverse mapping 

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

315 # arbitrary sky coordinates. 

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

317 if patchPolygon.intersection(detectorInPatchCoordinates): 

318 overlappingArea += patchPolygon.intersectionSingle( 

319 detectorInPatchCoordinates 

320 ).calculateArea() 

321 self.log.info( 

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

323 dataId["tract"], 

324 dataId["patch"], 

325 ) 

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

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

328 

329 if not overlappingArea: 

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

331 

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

333 

334 @timeMethod 

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

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

337 subtract from a science image. 

338 

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

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

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

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

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

344 CoaddPsf evaluated at the center of the detector. 

345 

346 Parameters 

347 ---------- 

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

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

350 `lsst.afw.image.Exposure`]] 

351 Coadds to be mosaicked, indexed on tract id. 

352 bbox : `lsst.geom.Box2I` 

353 Template Bounding box of the detector geometry onto which to 

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

355 template border. 

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

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

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

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

360 tract id. 

361 physical_filter : `str` 

362 Physical filter of the science image. 

363 visit : `int`, optional 

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

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

366 science image for the visit. 

367 

368 Returns 

369 ------- 

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

371 A struct with attributes: 

372 

373 ``template`` 

374 A template coadd exposure assembled out of patches, with 

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

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

377 

378 Raises 

379 ------ 

380 NoWorkFound 

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

382 """ 

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

384 

385 bbox.grow(self.config.templateBorderSize) 

386 

387 warped = {} 

388 catalogs = [] 

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

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

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

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

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

394 # different instruments are not supported. 

395 plateScaleFactor = None 

396 if self.config.doCorrectVariancePlateScale and coaddExposureHandles: 

397 plateScaleFactor = self._plateScaleFactor( 

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

399 ) 

400 

401 for tract in coaddExposureHandles: 

402 maskedImages, catalog, totalBox = self._makeExposureCatalog( 

403 coaddExposureHandles[tract], dataIds[tract] 

404 ) 

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

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

407 # Combine images from individual patches together. 

408 unwarped, count, included = self._merge( 

409 maskedImages, warpedBox, catalog[0].wcs 

410 ) 

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

412 del maskedImages 

413 if count == 0: 

414 self.log.info( 

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

416 tract, 

417 ) 

418 continue 

419 warpedBox.clip(totalBox) 

420 potentialInput = self.warper.warpExposure( 

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

422 ) 

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

424 del unwarped 

425 if np.all( 

426 potentialInput.mask.array 

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

428 ): 

429 self.log.info( 

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

431 tract, 

432 ) 

433 continue 

434 

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

436 tempCatalog = afwTable.ExposureCatalog(self.schema) 

437 tempCatalog.reserve(len(included)) 

438 for i in included: 

439 tempCatalog.append(catalog[i]) 

440 catalogs.append(tempCatalog) 

441 warped[tract] = potentialInput.maskedImage 

442 

443 if len(warped) == 0: 

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

445 

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

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

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

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

450 

451 if self.config.doScaleVariance: 

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

453 # needed. Note that the science variance is scaled 

454 # independently in ``AlardLuptonSubtractTask``. 

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

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

457 self.metadata["scaleTemplateVarianceFactor"] = varianceFactor 

458 

459 self._correctVariance(template, plateScaleFactor) 

460 

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

462 catalog = afwTable.ExposureCatalog(self.schema) 

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

464 for c in catalogs: 

465 catalog.extend(c) 

466 

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

468 self.checkHighVariance(template) 

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

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

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

472 identity = afwImage.PhotoCalib(1.0) 

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

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

475 photoCalib = identity 

476 template.setPhotoCalib(photoCalib) 

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

478 template.setPsf(self._makeGaussianPsf(template)) 

479 

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

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

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

483 template.getInfo().setCoaddInputs(coaddInputs) 

484 return pipeBase.Struct(template=template) 

485 

486 def _correctVariance(self, template, plateScaleFactor): 

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

488 

489 Compute a multiplicative correction factor for the template variance 

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

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

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

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

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

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

496 image pixel noise level would suggest. 

497 

498 Parameters 

499 ---------- 

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

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

502 plateScaleFactor : `float` or `None` 

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

504 

505 Raises 

506 ------ 

507 RuntimeError 

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

509 could not be reconstructed from the coadd inputs. 

510 """ 

511 scale = 1.0 

512 if self.config.doCorrectVariancePlateScale: 

513 if plateScaleFactor is None: 

514 raise RuntimeError( 

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

516 ) 

517 scale *= plateScaleFactor 

518 self.metadata["variancePlateScaleFactor"] = plateScaleFactor 

519 self.log.info( 

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

521 " coadd inputs.", 

522 plateScaleFactor, 

523 ) 

524 

525 if scale == 1.0: 

526 return 

527 template.variance.array *= scale 

528 self.metadata["templateVarianceCorrectionFactor"] = scale 

529 self.log.info( 

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

531 scale, 

532 ) 

533 

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

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

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

537 

538 Parameters 

539 ---------- 

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

541 `lsst.afw.image.Exposure` 

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

543 component is read; the pixels are left alone. 

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

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

546 bbox : `lsst.geom.Box2I` 

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

548 evaluate the science image pixel scale. 

549 

550 Returns 

551 ------- 

552 factor : `float` or `None` 

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

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

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

556 """ 

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

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

559 return None 

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

561 

562 for record in coaddInputs.ccds: 

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

564 # that gives a usable pixel area. 

565 originalWcs = record.getWcs() 

566 if originalWcs is None: 

567 continue 

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

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

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

571 factor = (pviScale/scienceScale)**2 

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

573 return factor 

574 return None 

575 

576 def checkHighVariance(self, template): 

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

578 

579 Parameters 

580 ---------- 

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

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

583 """ 

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

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

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

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

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

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

590 " pixels were unmasked for background estimation, but" 

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

592 else: 

593 varianceExposure = template.clone() 

594 varianceExposure.image.array = varianceExposure.variance.array 

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

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

597 highVariancePix = varianceBackground > threshold 

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

599 

600 @staticmethod 

601 def _checkInputs(dataIds, coaddExposures): 

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

603 the exposures all have the same photometric calibration. 

604 

605 Parameters 

606 ---------- 

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

608 Record of the tract and patch of each coaddExposure. 

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

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

611 `lsst.afw.image.Exposure` or 

612 `lsst.afw.image.Exposure`]] 

613 Coadds to be mosaicked. 

614 

615 Returns 

616 ------- 

617 band : `str` 

618 Filter band of all the input exposures. 

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

620 Photometric calibration of all of the input exposures. 

621 

622 Raises 

623 ------ 

624 RuntimeError 

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

626 all the same. 

627 """ 

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

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

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

631 band = bands.pop() 

632 photoCalibs = [ 

633 exposure.get(component="photoCalib") 

634 for exposures in coaddExposures.values() 

635 for exposure in exposures 

636 ] 

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

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

639 raise RuntimeError(msg) 

640 photoCalib = photoCalibs[0] 

641 return band, photoCalib 

642 

643 def _makeExposureCatalog(self, exposureRefs, dataIds): 

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

645 

646 Parameters 

647 ---------- 

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

649 `lsst.afw.image.Exposure`] 

650 Exposures to include in the catalog. 

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

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

653 "patch" entries. 

654 

655 Returns 

656 ------- 

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

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

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

660 Catalog of metadata for each exposure 

661 totalBox : `lsst.geom.Box2I` 

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

663 """ 

664 catalog = afwTable.ExposureCatalog(self.schema) 

665 catalog.reserve(len(exposureRefs)) 

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

667 images = {} 

668 totalBox = geom.Box2I() 

669 

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

671 images[dataId] = coadd.maskedImage 

672 bbox = coadd.getBBox() 

673 totalBox = totalBox.expandedTo(bbox) 

674 record = catalog.addNew() 

675 record.setPsf(coadd.psf) 

676 record.setWcs(coadd.wcs) 

677 record.setPhotoCalib(coadd.photoCalib) 

678 record.setBBox(bbox) 

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

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

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

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

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

684 record.set("weight", 1) 

685 

686 return images, catalog, totalBox 

687 

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

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

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

691 exposures. 

692 

693 Parameters 

694 ---------- 

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

696 `lsst.afw.image.Exposure`] 

697 Images to be merged into one larger bounding box. 

698 bbox : `lsst.geom.Box2I` 

699 Bounding box defining the image to merge into. 

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

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

702 

703 Returns 

704 ------- 

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

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

707 positions. 

708 count : `int` 

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

710 in the merged image. 

711 included : `list` [`int`] 

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

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

714 """ 

715 merged = afwImage.ExposureF(bbox, wcs) 

716 weights = afwImage.ImageF(bbox) 

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

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

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

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

721 clippedBox.clip(bbox) 

722 if clippedBox.area == 0: 

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

724 continue # nothing in this image overlaps the output 

725 maskedImage = maskedImage.subset(clippedBox) 

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

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

728 np.isfinite(maskedImage.variance.array) 

729 ) 

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

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

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

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

734 maskedImage.image.array[bad] = 0.0 

735 maskedImage.variance.array[bad] = 0.0 

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

737 maskedImage.mask.array[bad] = 0 

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

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

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

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

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

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

744 # Free memory before creating new large arrays 

745 del weight 

746 merged.maskedImage[clippedBox] += maskedImage 

747 included.append(i) 

748 

749 good = weights.array > 0 

750 

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

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

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

754 weights = weights.array[good] 

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

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

757 

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

759 

760 return merged, good.sum(), included 

761 

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

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

764 

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

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

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

768 went in to those pixels. 

769 

770 Parameters 

771 ---------- 

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

773 Generated template the PSF is for. 

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

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

776 of the input PSFs. 

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

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

779 

780 Returns 

781 ------- 

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

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

784 """ 

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

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

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

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

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

790 

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

792 coaddPsf = CoaddPsf( 

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

794 ) 

795 return coaddPsf 

796 

797 def _makeGaussianPsf(self, template): 

798 """Return a Gaussian PSF with the width of the template CoaddPsf. 

799 

800 Parameters 

801 ---------- 

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

803 Template with the CoaddPsf to approximate. The width is measured 

804 at the average position of the PSF, or averaged over a grid of 

805 positions if the PSF cannot be evaluated there. 

806 

807 Returns 

808 ------- 

809 psf : `lsst.afw.detection.GaussianPsf` 

810 Gaussian with the width of the template PSF. 

811 

812 Raises 

813 ------ 

814 ValueError 

815 Raised if the PSF cannot be evaluated at any grid position. 

816 """ 

817 # The CoaddPsf is not defined where the template has no inputs, which 

818 # can include its average position. 

819 try: 

820 fwhm = getPsfFwhm(template.psf) 

821 except lsst.pex.exceptions.Exception: 

822 # Catch a broad range of exceptions, since some are C++ only 

823 # Catching: 

824 # - lsst::geom::SingularTransformException 

825 # - lsst.pex.exceptions.InvalidParameterError 

826 # - lsst.pex.exceptions.RangeError 

827 self.log.info("Unable to evaluate PSF at the average position. " 

828 "Evaluating PSF on a grid of points.") 

829 fwhm = evaluateMeanPsfFwhm(template, 

830 fwhmExposureBuffer=self.config.fwhmExposureBuffer, 

831 fwhmExposureGrid=self.config.fwhmExposureGrid) 

832 sigma = fwhm/gaussian_sigma_to_fwhm 

833 size = 2*int(np.ceil(_PSF_STAMP_SIGMAS*sigma)) + 1 

834 self.metadata["templatePsfSigma"] = float(sigma) 

835 self.log.info("Approximated the template PSF with a Gaussian of sigma %.3f pixels.", sigma) 

836 return afwDetection.GaussianPsf(size, size, sigma) 

837 

838 

839class GetDcrTemplateConnections( 

840 GetTemplateConnections, 

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

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

843): 

844 visitInfo = pipeBase.connectionTypes.Input( 

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

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

847 storageClass="VisitInfo", 

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

849 ) 

850 dcrCoadds = pipeBase.connectionTypes.Input( 

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

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

853 storageClass="ExposureF", 

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

855 multiple=True, 

856 deferLoad=True, 

857 ) 

858 

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

860 super().__init__(config=config) 

861 self.inputs.remove("coaddExposures") 

862 

863 

864class GetDcrTemplateConfig( 

865 GetTemplateConfig, pipelineConnections=GetDcrTemplateConnections 

866): 

867 numSubfilters = pexConfig.Field( 

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

869 dtype=int, 

870 default=3, 

871 ) 

872 effectiveWavelength = pexConfig.Field( 

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

874 optional=False, 

875 dtype=float, 

876 ) 

877 bandwidth = pexConfig.Field( 

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

879 optional=False, 

880 dtype=float, 

881 ) 

882 

883 def validate(self): 

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

885 raise ValueError( 

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

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

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

889 ) 

890 

891 

892class GetDcrTemplateTask(GetTemplateTask): 

893 ConfigClass = GetDcrTemplateConfig