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

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 

23import functools 

24 

25import astropy.units 

26from astropy.stats import gaussian_sigma_to_fwhm 

27import numpy as np 

28 

29import lsst.afw.detection as afwDetection 

30import lsst.afw.image as afwImage 

31import lsst.geom as geom 

32import lsst.afw.geom as afwGeom 

33from lsst.afw.image import VisitInfo 

34import lsst.afw.table as afwTable 

35from lsst.afw.math._warper import computeWarpedBBox 

36import lsst.afw.math as afwMath 

37import lsst.pex.config as pexConfig 

38import lsst.pex.exceptions 

39import lsst.pipe.base as pipeBase 

40 

41from lsst.images import ( 

42 DifferenceImage, 

43 DifferenceImageTemplateInfo, 

44 get_legacy_optional_mask_planes, 

45 get_legacy_template_mask_planes, 

46) 

47from lsst.skymap import BaseSkyMap 

48from lsst.ip.diffim.dcrModel import DcrModel 

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

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

51from lsst.utils.timer import timeMethod 

52 

53__all__ = [ 

54 "GetTemplateTask", 

55 "GetTemplateConfig", 

56 "GetDcrTemplateTask", 

57 "GetDcrTemplateConfig", 

58] 

59 

60_PSF_STAMP_SIGMAS = 5 

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

62 

63_UNKNOWN_INSTRUMENT = "UNKNOWN" 

64"""Instrument recorded for a coadd or template made from observations of an 

65unknown instrument. 

66""" 

67 

68 

69@functools.cache 

70def _warpKernelSumSq(kernelName, nPhase=100): 

71 """Compute the phase-averaged sum of squared values of a warping kernel. 

72 

73 Warping with a kernel normalized to unit sum multiplies the variance by 

74 this sum, which is less than one. The power the warp took off the 

75 diagonal became covariance between neighboring pixels, which a variance 

76 plane cannot represent, so dividing by it restores the noise power: the 

77 variance of a flux measured over a region and not a single pixel. 

78 

79 Parameters 

80 ---------- 

81 kernelName : `str` 

82 Name of the warping kernel, as `lsst.afw.math.Warper` takes it. 

83 nPhase : `int`, optional 

84 Number of sub-pixel phases to average over. 

85 

86 Returns 

87 ------- 

88 sumSq : `float` 

89 Sum of the squared 2-D kernel values, averaged over sub-pixel 

90 phase. Converges with large ``nPhase`` to 0.798 for ``lanczos3``, 

91 4/9 for ``bilinear``, and 1 for no interpolation. 

92 

93 Notes 

94 ----- 

95 The kernels are separable and the phases of the two axes independent, 

96 so this is the square of the 1-D average. The warper normalizes the 

97 kernel to unit sum at each phase, which the analytic Lanczos integral 

98 does not include, so the average is computed from the kernel itself. 

99 """ 

100 # Disable the cache since we are evaluating fewer phase points than the 

101 # default cache size. 

102 kernel = afwMath.Warper(kernelName, cacheSize=0).getWarpingKernel() 

103 image = afwImage.ImageD(kernel.getDimensions()) 

104 total = 0.0 

105 for phase in (np.arange(nPhase) + 0.5)/nPhase: 

106 kernel.setKernelParameters((float(phase), 0.0)) 

107 kernel.computeImage(image, doNormalize=True) 

108 total += float((image.array.sum(axis=0)**2).sum()) 

109 return (total/nPhase)**2 

110 

111 

112class GetTemplateConnections( 

113 pipeBase.PipelineTaskConnections, 

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

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

116): 

117 bbox = pipeBase.connectionTypes.Input( 

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

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

120 storageClass="Box2I", 

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

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

123 ) 

124 detector = pipeBase.connectionTypes.Input( 

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

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

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

128 storageClass="Detector", 

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

130 ) 

131 obs_info = pipeBase.connectionTypes.Input( 

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

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

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

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

136 storageClass="ObservationInfo", 

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

138 ) 

139 wcs = pipeBase.connectionTypes.Input( 

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

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

142 storageClass="Wcs", 

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

144 ) 

145 skyMap = pipeBase.connectionTypes.Input( 

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

147 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

148 dimensions=("skymap",), 

149 storageClass="SkyMap", 

150 ) 

151 coaddExposures = pipeBase.connectionTypes.Input( 

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

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

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

155 storageClass="ExposureF", 

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

157 multiple=True, 

158 deferLoad=True, 

159 deferGraphConstraint=True, 

160 ) 

161 

162 template = pipeBase.connectionTypes.Output( 

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

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

165 storageClass="ExposureF", 

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

167 ) 

168 

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

170 super().__init__(config=config) 

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

172 del self.bbox 

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

174 self.coaddExposures = dataclasses.replace( 

175 self.coaddExposures, 

176 deferGraphConstraint=False, 

177 ) 

178 if config.output_image_type == "future": 

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

180 else: 

181 del self.obs_info 

182 

183 

184class GetTemplateConfig( 

185 pipeBase.PipelineTaskConfig, pipelineConnections=GetTemplateConnections 

186): 

187 templateBorderSize = pexConfig.Field( 

188 dtype=int, 

189 default=20, 

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

191 ) 

192 warp = pexConfig.ConfigField( 

193 dtype=afwMath.Warper.ConfigClass, 

194 doc="warper configuration", 

195 ) 

196 doCorrectVarianceWarpCorrelation = pexConfig.Field( 

197 dtype=bool, 

198 default=True, 

199 doc="Correct the template variance for the correlation that warping" 

200 " introduced between neighboring pixels? This covers both warps: the" 

201 " one that built the coadds, named by coaddWarpKernel, and the one" 

202 " this task performs, named by warp.warpingKernelName. It is a single" 

203 " scalar factor, fixed entirely by those two kernels and independent" 

204 " of the pixel scales involved.", 

205 ) 

206 coaddWarpKernel = pexConfig.ChoiceField( 

207 dtype=str, 

208 default="lanczos3", 

209 allowed={ 

210 "bilinear": "bilinear interpolation", 

211 "lanczos3": "Lanczos kernel of order 3", 

212 "lanczos4": "Lanczos kernel of order 4", 

213 "lanczos5": "Lanczos kernel of order 5", 

214 }, 

215 doc="Warping kernel that was used to build the input coadds. This" 

216 " cannot be recovered from the coadds themselves, and is used only if" 

217 " doCorrectVarianceWarpCorrelation is True.", 

218 ) 

219 coaddPsf = pexConfig.ConfigField( 

220 doc="Configuration for CoaddPsf", 

221 dtype=CoaddPsfConfig, 

222 ) 

223 fwhmExposureGrid = pexConfig.Field( 

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

225 dtype=int, 

226 default=10, 

227 ) 

228 fwhmExposureBuffer = pexConfig.Field( 

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

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

231 dtype=float, 

232 default=0.05, 

233 ) 

234 varianceBackground = pexConfig.ConfigurableField( 

235 target=SubtractBackgroundTask, 

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

237 ) 

238 highVarianceThreshold = pexConfig.RangeField( 

239 dtype=float, 

240 default=4, 

241 min=1, 

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

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

244 ) 

245 highVarianceMaskFraction = pexConfig.Field( 

246 dtype=float, 

247 default=0.1, 

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

249 " HIGH_VARIANCE mask plane.", 

250 ) 

251 doScaleVariance = pexConfig.Field( 

252 dtype=bool, 

253 default=False, 

254 doc="Scale the variance of the template image to match its per-pixel" 

255 " scatter? Warping correlates the noise, so when combined with" 

256 " doCorrectVarianceWarpCorrelation this overcounts the correction.", 

257 ) 

258 scaleVariance = pexConfig.ConfigurableField( 

259 target=ScaleVarianceTask, 

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

261 ) 

262 requireCoaddAtGraphBuild = pexConfig.Field( 

263 dtype=bool, 

264 default=False, 

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

266 " initial butler query during QuantumGraph generation.", 

267 ) 

268 output_image_type = pexConfig.ChoiceField[str]( 

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

270 allowed={ 

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

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

273 }, 

274 optional=False, 

275 default="legacy", 

276 ) 

277 raiseOnUndefinedMaskMap = pexConfig.Field( 

278 dtype=bool, 

279 default=False, 

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

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

282 " plane from the output instead." 

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

284 ) 

285 

286 def setDefaults(self): 

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

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

289 self.warp.cacheSize = 100000 

290 self.coaddPsf.cacheSize = self.warp.cacheSize 

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

292 # interpolation length for WCS evaluations. 

293 self.warp.interpLength = 100 

294 self.warp.warpingKernelName = "lanczos3" 

295 self.coaddPsf.warpingKernelName = self.warp.warpingKernelName 

296 

297 # Background subtraction of the variance plane 

298 self.varianceBackground.algorithm = "LINEAR" 

299 self.varianceBackground.binSize = 32 

300 self.varianceBackground.useApprox = False 

301 self.varianceBackground.statisticsProperty = "MEDIAN" 

302 self.varianceBackground.doFilterSuperPixels = True 

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

304 "EDGE", 

305 "DETECTED", 

306 "DETECTED_NEGATIVE", 

307 "NO_DATA", 

308 ] 

309 

310 

311class GetTemplateTask(pipeBase.PipelineTask): 

312 ConfigClass = GetTemplateConfig 

313 _DefaultName = "getTemplate" 

314 

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

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

317 if self.config.doScaleVariance: 

318 self.makeSubtask("scaleVariance") 

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

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

321 self.schema.addField( 

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

323 ) 

324 self.schema.addField( 

325 "patch", 

326 type=np.int32, 

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

328 ) 

329 self.schema.addField( 

330 "weight", 

331 type=float, 

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

333 ) 

334 self.schema.addField( 

335 "instrument", 

336 type=str, 

337 size=64, 

338 doc="Instrument whose observations were coadded to make the template, or UNKNOWN.", 

339 ) 

340 self.makeSubtask("varianceBackground") 

341 

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

343 inputs = butlerQC.get(inputRefs) 

344 detector = inputs.pop("detector") 

345 bbox = detector.getBBox() 

346 wcs = inputs.pop("wcs") 

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

348 coaddExposures = inputs.pop("coaddExposures") 

349 skymap = inputs.pop("skyMap") 

350 

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

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

353 

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

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

356 outputs = self.run( 

357 coaddExposureHandles=results.coaddExposures, 

358 bbox=bbox, 

359 wcs=wcs, 

360 dataIds=results.dataIds, 

361 physical_filter=physical_filter, 

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

363 instrument=butlerQC.quantum.dataId["instrument"], 

364 ) 

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

366 dataId = butlerQC.quantum.dataId 

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

368 dataId.universe) 

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

370 inputRefs.coaddExposures, 

371 detector=detector, 

372 exposureRecord=exposureRecord) 

373 butlerQC.put(outputs, outputRefs) 

374 

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

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

377 

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

379 instance and records the coadds contributed to it, along with the 

380 instrument they were made from. The Gaussian PSF set by `run` is 

381 converted to an `lsst.images.psfs.GaussianPointSpreadFunction`. 

382 

383 Parameters 

384 ---------- 

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

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

387 and its mask is modified in place. 

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

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

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

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

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

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

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

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

396 for, which supplies the observation metadata recorded on the 

397 output. 

398 

399 Returns 

400 ------- 

401 convertedTemplate : `lsst.images.DifferenceImage` 

402 The converted image. 

403 

404 Raises 

405 ------ 

406 RuntimeError 

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

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

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

410 spans more than one skymap. 

411 """ 

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

413 # template inputs to record which coadds contributed. 

414 catalog = template.getInfo().getCoaddInputs().ccds 

415 coaddPsf = self._makePsf(template, catalog, template.getWcs()) 

416 # `run` records the template's instrument on every input. 

417 instrument = catalog[0]["instrument"] 

418 template.setDetector(detector) 

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

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

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

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

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

424 for ref in coadd_refs} 

425 planeMap = get_legacy_template_mask_planes() 

426 if not self.config.raiseOnUndefinedMaskMap: 

427 self._clearUnmappedMaskPlanes(template.mask, planeMap) 

428 convertedTemplate = DifferenceImage.from_legacy( 

429 template, 

430 exposure_record=exposureRecord, 

431 unit=astropy.units.nJy, 

432 plane_map=planeMap, 

433 ) 

434 del template 

435 convertedTemplate.templates = DifferenceImageTemplateInfo.from_legacy_psf( 

436 convertedTemplate.sky_projection.pixel_frame, 

437 coaddPsf, 

438 skymaps.pop(), 

439 butler_info, 

440 log=self.log, 

441 instrument=instrument, 

442 ) 

443 return convertedTemplate 

444 

445 def _clearUnmappedMaskPlanes(self, mask, planeMap): 

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

447 

448 Parameters 

449 ---------- 

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

451 Template mask to modify in place. 

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

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

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

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

456 """ 

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

458 for name in mask.getMaskPlaneDict(): 

459 if name in keep: 

460 continue 

461 bitmask = mask.getPlaneBitMask(name) 

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

463 if nSet: 

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

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

466 mask.array &= ~bitmask 

467 

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

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

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

471 structure of their dataIds. 

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

473 

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

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

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

477 all possibly-matching coadd exposures. 

478 

479 Parameters 

480 ---------- 

481 coaddExposureHandles : `iterable` \ 

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

483 `lsst.afw.image.Exposure`] 

484 Dataset handles to exposures that might overlap the desired 

485 region. 

486 bbox : `lsst.geom.Box2I` 

487 Template bounding box of the pixel geometry onto which the 

488 coaddExposures will be resampled. 

489 skymap : `lsst.skymap.SkyMap` 

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

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

492 Template WCS onto which the coadds will be resampled. 

493 

494 Returns 

495 ------- 

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

497 A struct with attributes: 

498 

499 ``coaddExposures`` 

500 Dict of coadd exposures that overlap the projected bbox, 

501 indexed on tract id 

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

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

504 ``dataIds`` 

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

506 projected bbox, indexed on tract id 

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

508 

509 Raises 

510 ------ 

511 NoWorkFound 

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

513 WCS is None. 

514 """ 

515 if wcs is None: 

516 raise pipeBase.NoWorkFound( 

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

518 ) 

519 

520 # Exposure's validPolygon would be more accurate 

521 detectorPolygon = geom.Box2D(bbox) 

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

523 overlappingArea = 0 

524 coaddExposures = collections.defaultdict(list) 

525 dataIds = collections.defaultdict(list) 

526 

527 for coaddRef in coaddExposureHandles: 

528 dataId = coaddRef.dataId 

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

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

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

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

533 # detector coordinates because the skymap's inverse mapping 

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

535 # arbitrary sky coordinates. 

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

537 if patchPolygon.intersection(detectorInPatchCoordinates): 

538 overlappingArea += patchPolygon.intersectionSingle( 

539 detectorInPatchCoordinates 

540 ).calculateArea() 

541 self.log.info( 

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

543 dataId["tract"], 

544 dataId["patch"], 

545 ) 

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

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

548 

549 if not overlappingArea: 

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

551 

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

553 

554 @timeMethod 

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

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

557 subtract from a science image. 

558 

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

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

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

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

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

564 CoaddPsf evaluated at the center of the detector. 

565 

566 The template is built from coadds of one instrument: the science 

567 image's ``instrument`` if any coadd was made from it, and otherwise the 

568 one that most of the coadds were made from. Coadds made from another 

569 instrument are dropped, and so are coadds that record no instrument 

570 unless none of the coadds records one. 

571 

572 Parameters 

573 ---------- 

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

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

576 `lsst.afw.image.Exposure`]] 

577 Coadds to be mosaicked, indexed on tract id. 

578 bbox : `lsst.geom.Box2I` 

579 Input Bounding box of the detector geometry onto which to 

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

581 will be grown by ``templateBorderSize``. 

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

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

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

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

586 tract id. 

587 physical_filter : `str` 

588 Physical filter of the science image. 

589 visit : `int`, optional 

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

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

592 science image for the visit. 

593 instrument : `str`, optional 

594 Instrument of the science image. Coadds made from it are used in 

595 preference to coadds from any other instrument. 

596 

597 Returns 

598 ------- 

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

600 A struct with attributes: 

601 

602 ``template`` 

603 A template coadd exposure assembled out of patches, with 

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

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

606 

607 Raises 

608 ------ 

609 NoWorkFound 

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

611 RuntimeError 

612 If the template is not known to be from the science image's 

613 instrument and no coadd records the plate scale of its input 

614 images, so the template variance cannot be corrected. 

615 """ 

616 band, photoCalib, templateInstrument = self._checkInputs(dataIds, coaddExposureHandles, instrument) 

617 

618 bbox = geom.Box2I(bbox) 

619 bbox.grow(self.config.templateBorderSize) 

620 

621 warped = {} 

622 catalogs = [] 

623 plateScaleFactor = 1.0 if instrument == templateInstrument else None 

624 

625 for tract in coaddExposureHandles: 

626 maskedImages, catalog, totalBox, coaddInputs = self._makeExposureCatalog( 

627 coaddExposureHandles[tract], dataIds[tract], templateInstrument 

628 ) 

629 if len(catalog) == 0: 

630 # Every coadd in this tract was dropped for its instrument. 

631 continue 

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

633 # of the science image. This will only be different from 1 if the 

634 # coadd was from a different instrument. A single correction factor 

635 # is calculated from the first usable coadd and applied to all. 

636 if plateScaleFactor is None: 

637 plateScaleFactor = self._plateScaleFactor(coaddInputs, wcs, bbox) 

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

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

640 # Combine images from individual patches together. 

641 unwarped, count, included = self._merge( 

642 maskedImages, warpedBox, catalog[0].wcs 

643 ) 

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

645 del maskedImages 

646 if count == 0: 

647 self.log.info( 

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

649 tract, 

650 ) 

651 continue 

652 warpedBox.clip(totalBox) 

653 potentialInput = self.warper.warpExposure( 

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

655 ) 

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

657 del unwarped 

658 if np.all( 

659 potentialInput.mask.array 

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

661 ): 

662 self.log.info( 

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

664 tract, 

665 ) 

666 continue 

667 

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

669 tempCatalog = afwTable.ExposureCatalog(self.schema) 

670 tempCatalog.reserve(len(included)) 

671 for i in included: 

672 tempCatalog.append(catalog[i]) 

673 catalogs.append(tempCatalog) 

674 warped[tract] = potentialInput.maskedImage 

675 

676 if len(warped) == 0: 

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

678 

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

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

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

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

683 

684 if self.config.doScaleVariance: 

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

686 # needed. Note that the science variance is scaled 

687 # independently in ``AlardLuptonSubtractTask``. 

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

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

690 self.metadata["scaleTemplateVarianceFactor"] = varianceFactor 

691 

692 self._correctVariance(template, plateScaleFactor) 

693 

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

695 catalog = afwTable.ExposureCatalog(self.schema) 

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

697 for c in catalogs: 

698 catalog.extend(c) 

699 

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

701 self.checkHighVariance(template) 

702 if visit is not None: 

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

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

705 identity = afwImage.PhotoCalib(1.0) 

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

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

708 photoCalib = identity 

709 template.setPhotoCalib(photoCalib) 

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

711 template.setPsf(self._makeGaussianPsf(template)) 

712 

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

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

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

716 template.getInfo().setCoaddInputs(coaddInputs) 

717 return pipeBase.Struct(template=template) 

718 

719 def _correctVariance(self, template, plateScaleFactor): 

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

721 

722 Warping introduces a covariance between neighboring pixels, while the 

723 input variance plane captures only the diagonal term. Compute a 

724 multiplicative correction factor to apply to the final variance so that 

725 the detection significance on the convolved difference image is correct 

726 on average. 

727 

728 Separately, also apply a scalar multiplicative correction factor to the 

729 template variance if the plate scale of the coadd's constituent images 

730 is different than the science image the template is being constructed 

731 for. This should only be necessary if the coadd images were from a 

732 different instrument than the science image, since the scale change 

733 otherwise cancels between the two warps. If the plate scale of the 

734 science instrument is smaller than the plate scale of the coadd 

735 instrument, then the template pixels will be correlated and the true 

736 variance will be higher than the image pixel noise level would suggest. 

737 

738 Parameters 

739 ---------- 

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

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

742 plateScaleFactor : `float` or `None` 

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

744 

745 Raises 

746 ------ 

747 RuntimeError 

748 Raised if ``plateScaleFactor`` is `None`. 

749 """ 

750 if plateScaleFactor is None: 

751 raise RuntimeError( 

752 "Could not determine the plate scale of the images the template coadds were made from," 

753 " so the template variance cannot be corrected for the change in pixel area." 

754 ) 

755 self.metadata["variancePlateScaleFactor"] = plateScaleFactor 

756 correlationFactor = 1.0 

757 if self.config.doCorrectVarianceWarpCorrelation: 

758 # Both warps correlate the noise: the one that built the coadds, 

759 # and the one this task just performed. 

760 correlationFactor = 1/(_warpKernelSumSq(self.config.coaddWarpKernel) 

761 * _warpKernelSumSq(self.config.warp.warpingKernelName)) 

762 self.metadata["varianceWarpCorrelationFactor"] = correlationFactor 

763 scale = correlationFactor*plateScaleFactor 

764 if scale == 1.0: 

765 return 

766 template.variance.array *= scale 

767 self.metadata["templateVarianceCorrectionFactor"] = scale 

768 self.log.info( 

769 "Variance correction factor due to plate scale changes: %.4f", 

770 plateScaleFactor, 

771 ) 

772 self.log.info( 

773 "Variance correction factor due to warping correlations: %.4f", 

774 correlationFactor, 

775 ) 

776 self.log.info( 

777 "Overall template variance correction factor: %.4f", 

778 scale, 

779 ) 

780 

781 def _plateScaleFactor(self, coaddInputs, wcs, bbox): 

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

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

784 

785 Parameters 

786 ---------- 

787 coaddInputs : `lsst.afw.image.CoaddInputs` or `None` 

788 Input records of one of the coadd patches. 

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

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

791 bbox : `lsst.geom.Box2I` 

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

793 evaluate the science image pixel scale. 

794 

795 Returns 

796 ------- 

797 factor : `float` or `None` 

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

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

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

801 """ 

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

803 return None 

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

805 

806 for record in coaddInputs.ccds: 

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

808 # that gives a usable pixel area. 

809 originalWcs = record.getWcs() 

810 if originalWcs is None: 

811 continue 

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

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

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

815 factor = (pviScale/scienceScale)**2 

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

817 return factor 

818 return None 

819 

820 def checkHighVariance(self, template): 

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

822 

823 Parameters 

824 ---------- 

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

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

827 """ 

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

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

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

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

832 if goodFraction < self.config.highVarianceMaskFraction: 

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

834 " pixels were unmasked for background estimation, but" 

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

836 else: 

837 varianceExposure = template.clone() 

838 varianceExposure.image.array = varianceExposure.variance.array 

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

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

841 highVariancePix = varianceBackground > threshold 

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

843 

844 @staticmethod 

845 def _checkInputs(dataIds, coaddExposures, scienceInstrument=None): 

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

847 the exposures all have the same photometric calibration, and choose 

848 the instrument of the template. 

849 

850 Parameters 

851 ---------- 

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

853 Record of the tract and patch of each coaddExposure. 

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

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

856 `lsst.afw.image.Exposure` or 

857 `lsst.afw.image.Exposure`]] 

858 Coadds to be mosaicked. 

859 scienceInstrument : `str`, optional 

860 Instrument of the science image. 

861 

862 Returns 

863 ------- 

864 band : `str` 

865 Filter band of all the input exposures. 

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

867 Photometric calibration of all of the input exposures. 

868 instrument : `str` 

869 ``scienceInstrument`` if any coadd was made from it, otherwise the 

870 instrument that the most coadds were made from, or ``"UNKNOWN"`` if 

871 no coadd records one. A tie goes to the instrument whose name sorts 

872 first. 

873 

874 Raises 

875 ------ 

876 RuntimeError 

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

878 all the same. 

879 """ 

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

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

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

883 band = bands.pop() 

884 photoCalibs = [] 

885 instrumentCounts = collections.Counter() 

886 for exposures in coaddExposures.values(): 

887 for exposure in exposures: 

888 # Read both components of a coadd in turn, so a butler cache 

889 # can serve the second read from the first download. 

890 photoCalibs.append(exposure.get(component="photoCalib")) 

891 instrumentCounts[GetTemplateTask._getInstrument(exposure.get(component="coaddInputs"))] += 1 

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

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

894 raise RuntimeError(msg) 

895 photoCalib = photoCalibs[0] 

896 # Coadds that record no instrument do not count toward the choice. 

897 del instrumentCounts[_UNKNOWN_INSTRUMENT] 

898 if scienceInstrument in instrumentCounts: 

899 return band, photoCalib, scienceInstrument 

900 if not instrumentCounts: 

901 return band, photoCalib, _UNKNOWN_INSTRUMENT 

902 instrument = max(sorted(instrumentCounts), key=instrumentCounts.get)