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

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_UNKNOWN_INSTRUMENT = "UNKNOWN" 

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

64unknown instrument. 

65""" 

66 

67 

68class GetTemplateConnections( 

69 pipeBase.PipelineTaskConnections, 

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

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

72): 

73 bbox = pipeBase.connectionTypes.Input( 

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

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

76 storageClass="Box2I", 

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

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

79 ) 

80 detector = pipeBase.connectionTypes.Input( 

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

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

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

84 storageClass="Detector", 

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

86 ) 

87 obs_info = pipeBase.connectionTypes.Input( 

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

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

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

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

92 storageClass="ObservationInfo", 

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

94 ) 

95 wcs = pipeBase.connectionTypes.Input( 

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

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

98 storageClass="Wcs", 

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

100 ) 

101 skyMap = pipeBase.connectionTypes.Input( 

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

103 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

104 dimensions=("skymap",), 

105 storageClass="SkyMap", 

106 ) 

107 coaddExposures = pipeBase.connectionTypes.Input( 

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

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

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

111 storageClass="ExposureF", 

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

113 multiple=True, 

114 deferLoad=True, 

115 deferGraphConstraint=True, 

116 ) 

117 

118 template = pipeBase.connectionTypes.Output( 

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

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

121 storageClass="ExposureF", 

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

123 ) 

124 

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

126 super().__init__(config=config) 

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

128 del self.bbox 

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

130 self.coaddExposures = dataclasses.replace( 

131 self.coaddExposures, 

132 deferGraphConstraint=False, 

133 ) 

134 if config.output_image_type == "future": 

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

136 else: 

137 del self.obs_info 

138 

139 

140class GetTemplateConfig( 

141 pipeBase.PipelineTaskConfig, pipelineConnections=GetTemplateConnections 

142): 

143 templateBorderSize = pexConfig.Field( 

144 dtype=int, 

145 default=20, 

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

147 ) 

148 warp = pexConfig.ConfigField( 

149 dtype=afwMath.Warper.ConfigClass, 

150 doc="warper configuration", 

151 ) 

152 coaddPsf = pexConfig.ConfigField( 

153 doc="Configuration for CoaddPsf", 

154 dtype=CoaddPsfConfig, 

155 ) 

156 fwhmExposureGrid = pexConfig.Field( 

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

158 dtype=int, 

159 default=10, 

160 ) 

161 fwhmExposureBuffer = pexConfig.Field( 

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

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

164 dtype=float, 

165 default=0.05, 

166 ) 

167 varianceBackground = pexConfig.ConfigurableField( 

168 target=SubtractBackgroundTask, 

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

170 ) 

171 highVarianceThreshold = pexConfig.RangeField( 

172 dtype=float, 

173 default=4, 

174 min=1, 

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

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

177 ) 

178 highVarianceMaskFraction = pexConfig.Field( 

179 dtype=float, 

180 default=0.1, 

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

182 " HIGH_VARIANCE mask plane.", 

183 ) 

184 doScaleVariance = pexConfig.Field( 

185 dtype=bool, 

186 default=True, 

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

188 ) 

189 scaleVariance = pexConfig.ConfigurableField( 

190 target=ScaleVarianceTask, 

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

192 ) 

193 requireCoaddAtGraphBuild = pexConfig.Field( 

194 dtype=bool, 

195 default=False, 

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

197 " initial butler query during QuantumGraph generation.", 

198 ) 

199 output_image_type = pexConfig.ChoiceField[str]( 

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

201 allowed={ 

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

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

204 }, 

205 optional=False, 

206 default="legacy", 

207 ) 

208 raiseOnUndefinedMaskMap = pexConfig.Field( 

209 dtype=bool, 

210 default=False, 

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

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

213 " plane from the output instead." 

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

215 ) 

216 

217 def setDefaults(self): 

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

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

220 self.warp.cacheSize = 100000 

221 self.coaddPsf.cacheSize = self.warp.cacheSize 

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

223 # interpolation length for WCS evaluations. 

224 self.warp.interpLength = 100 

225 self.warp.warpingKernelName = "lanczos3" 

226 self.coaddPsf.warpingKernelName = self.warp.warpingKernelName 

227 

228 # Background subtraction of the variance plane 

229 self.varianceBackground.algorithm = "LINEAR" 

230 self.varianceBackground.binSize = 32 

231 self.varianceBackground.useApprox = False 

232 self.varianceBackground.statisticsProperty = "MEDIAN" 

233 self.varianceBackground.doFilterSuperPixels = True 

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

235 "EDGE", 

236 "DETECTED", 

237 "DETECTED_NEGATIVE", 

238 "NO_DATA", 

239 ] 

240 

241 

242class GetTemplateTask(pipeBase.PipelineTask): 

243 ConfigClass = GetTemplateConfig 

244 _DefaultName = "getTemplate" 

245 

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

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

248 if self.config.doScaleVariance: 

249 self.makeSubtask("scaleVariance") 

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

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

252 self.schema.addField( 

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

254 ) 

255 self.schema.addField( 

256 "patch", 

257 type=np.int32, 

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

259 ) 

260 self.schema.addField( 

261 "weight", 

262 type=float, 

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

264 ) 

265 self.schema.addField( 

266 "instrument", 

267 type=str, 

268 size=64, 

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

270 ) 

271 self.makeSubtask("varianceBackground") 

272 

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

274 inputs = butlerQC.get(inputRefs) 

275 detector = inputs.pop("detector") 

276 bbox = detector.getBBox() 

277 wcs = inputs.pop("wcs") 

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

279 coaddExposures = inputs.pop("coaddExposures") 

280 skymap = inputs.pop("skyMap") 

281 

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

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

284 

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

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

287 outputs = self.run( 

288 coaddExposureHandles=results.coaddExposures, 

289 bbox=bbox, 

290 wcs=wcs, 

291 dataIds=results.dataIds, 

292 physical_filter=physical_filter, 

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

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

295 ) 

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

297 dataId = butlerQC.quantum.dataId 

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

299 dataId.universe) 

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

301 inputRefs.coaddExposures, 

302 detector=detector, 

303 exposureRecord=exposureRecord) 

304 butlerQC.put(outputs, outputRefs) 

305 

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

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

308 

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

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

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

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

313 

314 Parameters 

315 ---------- 

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

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

318 and its mask is modified in place. 

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

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

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

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

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

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

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

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

327 for, which supplies the observation metadata recorded on the 

328 output. 

329 

330 Returns 

331 ------- 

332 convertedTemplate : `lsst.images.DifferenceImage` 

333 The converted image. 

334 

335 Raises 

336 ------ 

337 RuntimeError 

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

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

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

341 spans more than one skymap. 

342 """ 

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

344 # template inputs to record which coadds contributed. 

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

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

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

348 instrument = catalog[0]["instrument"] 

349 template.setDetector(detector) 

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

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

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

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

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

355 for ref in coadd_refs} 

356 planeMap = get_legacy_template_mask_planes() 

357 if not self.config.raiseOnUndefinedMaskMap: 

358 self._clearUnmappedMaskPlanes(template.mask, planeMap) 

359 convertedTemplate = DifferenceImage.from_legacy( 

360 template, 

361 exposure_record=exposureRecord, 

362 unit=astropy.units.nJy, 

363 plane_map=planeMap, 

364 ) 

365 del template 

366 convertedTemplate.templates = DifferenceImageTemplateInfo.from_legacy_psf( 

367 convertedTemplate.sky_projection.pixel_frame, 

368 coaddPsf, 

369 skymaps.pop(), 

370 butler_info, 

371 log=self.log, 

372 instrument=instrument, 

373 ) 

374 return convertedTemplate 

375 

376 def _clearUnmappedMaskPlanes(self, mask, planeMap): 

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

378 

379 Parameters 

380 ---------- 

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

382 Template mask to modify in place. 

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

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

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

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

387 """ 

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

389 for name in mask.getMaskPlaneDict(): 

390 if name in keep: 

391 continue 

392 bitmask = mask.getPlaneBitMask(name) 

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

394 if nSet: 

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

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

397 mask.array &= ~bitmask 

398 

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

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

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

402 structure of their dataIds. 

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

404 

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

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

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

408 all possibly-matching coadd exposures. 

409 

410 Parameters 

411 ---------- 

412 coaddExposureHandles : `iterable` \ 

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

414 `lsst.afw.image.Exposure`] 

415 Dataset handles to exposures that might overlap the desired 

416 region. 

417 bbox : `lsst.geom.Box2I` 

418 Template bounding box of the pixel geometry onto which the 

419 coaddExposures will be resampled. 

420 skymap : `lsst.skymap.SkyMap` 

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

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

423 Template WCS onto which the coadds will be resampled. 

424 

425 Returns 

426 ------- 

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

428 A struct with attributes: 

429 

430 ``coaddExposures`` 

431 Dict of coadd exposures that overlap the projected bbox, 

432 indexed on tract id 

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

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

435 ``dataIds`` 

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

437 projected bbox, indexed on tract id 

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

439 

440 Raises 

441 ------ 

442 NoWorkFound 

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

444 WCS is None. 

445 """ 

446 if wcs is None: 

447 raise pipeBase.NoWorkFound( 

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

449 ) 

450 

451 # Exposure's validPolygon would be more accurate 

452 detectorPolygon = geom.Box2D(bbox) 

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

454 overlappingArea = 0 

455 coaddExposures = collections.defaultdict(list) 

456 dataIds = collections.defaultdict(list) 

457 

458 for coaddRef in coaddExposureHandles: 

459 dataId = coaddRef.dataId 

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

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

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

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

464 # detector coordinates because the skymap's inverse mapping 

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

466 # arbitrary sky coordinates. 

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

468 if patchPolygon.intersection(detectorInPatchCoordinates): 

469 overlappingArea += patchPolygon.intersectionSingle( 

470 detectorInPatchCoordinates 

471 ).calculateArea() 

472 self.log.info( 

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

474 dataId["tract"], 

475 dataId["patch"], 

476 ) 

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

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

479 

480 if not overlappingArea: 

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

482 

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

484 

485 @timeMethod 

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

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

488 subtract from a science image. 

489 

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

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

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

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

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

495 CoaddPsf evaluated at the center of the detector. 

496 

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

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

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

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

501 unless none of the coadds records one. 

502 

503 Parameters 

504 ---------- 

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

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

507 `lsst.afw.image.Exposure`]] 

508 Coadds to be mosaicked, indexed on tract id. 

509 bbox : `lsst.geom.Box2I` 

510 Input Bounding box of the detector geometry onto which to 

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

512 will be grown by ``templateBorderSize``. 

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

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

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

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

517 tract id. 

518 physical_filter : `str` 

519 Physical filter of the science image. 

520 visit : `int`, optional 

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

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

523 science image for the visit. 

524 instrument : `str`, optional 

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

526 preference to coadds from any other instrument. 

527 

528 Returns 

529 ------- 

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

531 A struct with attributes: 

532 

533 ``template`` 

534 A template coadd exposure assembled out of patches, with 

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

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

537 

538 Raises 

539 ------ 

540 NoWorkFound 

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

542 RuntimeError 

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

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

545 images, so the template variance cannot be corrected. 

546 """ 

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

548 

549 bbox = geom.Box2I(bbox) 

550 bbox.grow(self.config.templateBorderSize) 

551 

552 warped = {} 

553 catalogs = [] 

554 plateScaleFactor = 1.0 if instrument == templateInstrument else None 

555 

556 for tract in coaddExposureHandles: 

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

558 coaddExposureHandles[tract], dataIds[tract], templateInstrument 

559 ) 

560 if len(catalog) == 0: 

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

562 continue 

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

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

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

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

567 if plateScaleFactor is None: 

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

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

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

571 # Combine images from individual patches together. 

572 unwarped, count, included = self._merge( 

573 maskedImages, warpedBox, catalog[0].wcs 

574 ) 

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

576 del maskedImages 

577 if count == 0: 

578 self.log.info( 

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

580 tract, 

581 ) 

582 continue 

583 warpedBox.clip(totalBox) 

584 potentialInput = self.warper.warpExposure( 

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

586 ) 

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

588 del unwarped 

589 if np.all( 

590 potentialInput.mask.array 

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

592 ): 

593 self.log.info( 

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

595 tract, 

596 ) 

597 continue 

598 

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

600 tempCatalog = afwTable.ExposureCatalog(self.schema) 

601 tempCatalog.reserve(len(included)) 

602 for i in included: 

603 tempCatalog.append(catalog[i]) 

604 catalogs.append(tempCatalog) 

605 warped[tract] = potentialInput.maskedImage 

606 

607 if len(warped) == 0: 

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

609 

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

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

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

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

614 

615 if self.config.doScaleVariance: 

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

617 # needed. Note that the science variance is scaled 

618 # independently in ``AlardLuptonSubtractTask``. 

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

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

621 self.metadata["scaleTemplateVarianceFactor"] = varianceFactor 

622 

623 self._correctVariance(template, plateScaleFactor) 

624 

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

626 catalog = afwTable.ExposureCatalog(self.schema) 

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

628 for c in catalogs: 

629 catalog.extend(c) 

630 

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

632 self.checkHighVariance(template) 

633 if visit is not None: 

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

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

636 identity = afwImage.PhotoCalib(1.0) 

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

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

639 photoCalib = identity 

640 template.setPhotoCalib(photoCalib) 

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

642 template.setPsf(self._makeGaussianPsf(template)) 

643 

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

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

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

647 template.getInfo().setCoaddInputs(coaddInputs) 

648 return pipeBase.Struct(template=template) 

649 

650 def _correctVariance(self, template, plateScaleFactor): 

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

652 

653 Apply a multiplicative correction factor to the template variance if 

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

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

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

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

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

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

660 image pixel noise level would suggest. 

661 

662 Parameters 

663 ---------- 

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

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

666 plateScaleFactor : `float` or `None` 

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

668 

669 Raises 

670 ------ 

671 RuntimeError 

672 Raised if ``plateScaleFactor`` is `None`. 

673 """ 

674 if plateScaleFactor is None: 

675 raise RuntimeError( 

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

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

678 ) 

679 self.metadata["variancePlateScaleFactor"] = plateScaleFactor 

680 if plateScaleFactor == 1.0: 

681 return 

682 template.variance.array *= plateScaleFactor 

683 self.metadata["templateVarianceCorrectionFactor"] = plateScaleFactor 

684 self.log.info( 

685 "Corrected the template variance plane for the change in pixel area by a factor of %.4f.", 

686 plateScaleFactor, 

687 ) 

688 

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

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

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

692 

693 Parameters 

694 ---------- 

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

696 Input records of one of the coadd patches. 

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

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

699 bbox : `lsst.geom.Box2I` 

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

701 evaluate the science image pixel scale. 

702 

703 Returns 

704 ------- 

705 factor : `float` or `None` 

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

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

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

709 """ 

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

711 return None 

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

713 

714 for record in coaddInputs.ccds: 

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

716 # that gives a usable pixel area. 

717 originalWcs = record.getWcs() 

718 if originalWcs is None: 

719 continue 

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

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

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

723 factor = (pviScale/scienceScale)**2 

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

725 return factor 

726 return None 

727 

728 def checkHighVariance(self, template): 

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

730 

731 Parameters 

732 ---------- 

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

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

735 """ 

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

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

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

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

740 if goodFraction < self.config.highVarianceMaskFraction: 

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

742 " pixels were unmasked for background estimation, but" 

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

744 else: 

745 varianceExposure = template.clone() 

746 varianceExposure.image.array = varianceExposure.variance.array 

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

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

749 highVariancePix = varianceBackground > threshold 

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

751 

752 @staticmethod 

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

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

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

756 the instrument of the template. 

757 

758 Parameters 

759 ---------- 

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

761 Record of the tract and patch of each coaddExposure. 

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

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

764 `lsst.afw.image.Exposure` or 

765 `lsst.afw.image.Exposure`]] 

766 Coadds to be mosaicked. 

767 scienceInstrument : `str`, optional 

768 Instrument of the science image. 

769 

770 Returns 

771 ------- 

772 band : `str` 

773 Filter band of all the input exposures. 

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

775 Photometric calibration of all of the input exposures. 

776 instrument : `str` 

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

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

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

780 first. 

781 

782 Raises 

783 ------ 

784 RuntimeError 

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

786 all the same. 

787 """ 

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

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

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

791 band = bands.pop() 

792 photoCalibs = [] 

793 instrumentCounts = collections.Counter() 

794 for exposures in coaddExposures.values(): 

795 for exposure in exposures: 

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

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

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

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

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

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

802 raise RuntimeError(msg) 

803 photoCalib = photoCalibs[0] 

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

805 del instrumentCounts[_UNKNOWN_INSTRUMENT] 

806 if scienceInstrument in instrumentCounts: 

807 return band, photoCalib, scienceInstrument 

808 if not instrumentCounts: 

809 return band, photoCalib, _UNKNOWN_INSTRUMENT 

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

811 return band, photoCalib, instrument 

812 

813 @staticmethod 

814 def _getInstrument(coaddInputs): 

815 """Return the instrument whose observations were coadded to make a 

816 coadd. 

817 

818 Parameters 

819 ---------- 

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

821 Input records of the coadd. 

822 

823 Returns 

824 ------- 

825 instrument : `str` 

826 Instrument label of the coadd, or ``"UNKNOWN"`` if it records none. 

827 

828 Notes 

829 ----- 

830 Every input of a coadd comes from the same instrument, so the first 

831 input record stands for all of them. 

832 """ 

833 if coaddInputs is None or len(coaddInputs.ccds) == 0: 

834 return _UNKNOWN_INSTRUMENT 

835 visitInfo = coaddInputs.ccds[0].getVisitInfo() 

836 label = visitInfo.getInstrumentLabel() if visitInfo is not None else "" 

837 return label or _UNKNOWN_INSTRUMENT 

838 

839 def _makeExposureCatalog(self, exposureRefs, dataIds, instrument=_UNKNOWN_INSTRUMENT): 

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

841 

842 Parameters 

843 ---------- 

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

845 `lsst.afw.image.Exposure`] 

846 Exposures to include in the catalog. 

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

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

849 "patch" entries. 

850 instrument : `str`, optional 

851 Instrument of the template, recorded in the catalog. Exposures 

852 made from another instrument are dropped, as are exposures that 

853 record no instrument unless ``instrument`` is ``"UNKNOWN"``. 

854 

855 Returns 

856 ------- 

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

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

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

860 Catalog of metadata for each exposure 

861 totalBox : `lsst.geom.Box2I` 

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

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

864 Input records of the first exposure kept, or `None` if none is 

865 kept. 

866 """ 

867 catalog = afwTable.ExposureCatalog(self.schema) 

868 catalog.reserve(len(exposureRefs)) 

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

870 images = {} 

871 totalBox = geom.Box2I() 

872 firstCoaddInputs = None 

873 

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

875 coaddInputs = coadd.getInfo().getCoaddInputs() 

876 coaddInstrument = self._getInstrument(coaddInputs) 

877 if coaddInstrument != instrument: 

878 self.log.info( 

879 "Dropping template coadd %s: it was made from %s, but the template is made from %s.", 

880 dataId, 

881 coaddInstrument, 

882 instrument, 

883 ) 

884 continue 

885 if not images: 

886 firstCoaddInputs = coaddInputs 

887 images[dataId] = coadd.maskedImage 

888 bbox = coadd.getBBox() 

889 totalBox = totalBox.expandedTo(bbox) 

890 record = catalog.addNew() 

891 record.setPsf(coadd.psf) 

892 record.setWcs(coadd.wcs) 

893 record.setPhotoCalib(coadd.photoCalib) 

894 record.setBBox(bbox) 

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

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