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

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 

37from lsst.daf.base import DateTime 

38import lsst.pex.config as pexConfig 

39import lsst.pex.exceptions 

40import lsst.pipe.base as pipeBase 

41 

42from lsst.images import ( 

43 Box, 

44 DifferenceImage, 

45 DifferenceImageTemplateInfo, 

46 DifferenceImageTemplateInput, 

47 GeneralFrame, 

48 NoOverlapError, 

49 Polygon, 

50 Transform, 

51 get_legacy_optional_mask_planes, 

52 get_legacy_template_mask_planes, 

53) 

54from lsst.skymap import BaseSkyMap 

55from lsst.ip.diffim.dcrModel import DcrModel 

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

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

58from lsst.utils.timer import timeMethod 

59 

60__all__ = [ 

61 "GetTemplateTask", 

62 "GetTemplateConfig", 

63 "GetDcrTemplateTask", 

64 "GetDcrTemplateConfig", 

65] 

66 

67_PSF_STAMP_SIGMAS = 5 

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

69 

70_UNKNOWN_INSTRUMENT = "UNKNOWN" 

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

72unknown instrument. 

73""" 

74 

75_FOOTPRINT_MARGIN = 32.0 

76"""Distance, in coadd pixels, by which the science image's footprint in a 

77coadd is grown before a coadd input is clipped to it. It must exceed twice 

78the distance by which a distorted science image edge departs from a straight 

79line between its corners, and limits how far outside the science image its 

80WCS is evaluated. 

81""" 

82 

83 

84def _pixelTransform(fromWcs, toWcs): 

85 """Return a chainable transformation between two pixel frames. 

86 

87 Parameters 

88 ---------- 

89 fromWcs, toWcs : `lsst.afw.geom.SkyWcs` 

90 WCSs of the input and output pixels. 

91 

92 Returns 

93 ------- 

94 transform : `lsst.images.Transform` 

95 Transform from ``fromWcs`` pixels to ``toWcs`` pixels. 

96 """ 

97 return Transform.from_legacy(afwGeom.makeWcsPairTransform(fromWcs, toWcs), 

98 GeneralFrame(unit=astropy.units.pix), 

99 GeneralFrame(unit=astropy.units.pix) 

100 ) 

101 

102 

103def _getImageOverlap(record, coaddWcs, scienceFootprint, coaddToScienceTransform, scienceBox, border): 

104 """Return the region enclosing the overlap of the science image and one 

105 input image of the coadd. 

106 

107 Note that there are three different WCS's involved here: the WCS of the 

108 original image that was warped to contribute to the coadd, the coadd WCS, 

109 and the current science image WCS that we are building a template for. 

110 

111 Parameters 

112 ---------- 

113 record : `lsst.afw.table.ExposureRecord` 

114 Input record of a coadd, with the input's WCS, bounding box and 

115 optional valid polygon. 

116 coaddWcs : `lsst.afw.geom.SkyWcs` 

117 WCS of the coadd. 

118 scienceFootprint : `lsst.images.Region` 

119 ``scienceBox`` in coadd pixels, grown by ``_FOOTPRINT_MARGIN``. 

120 coaddToScienceTransform : `lsst.images.Transform` 

121 Transform from coadd pixels to science image pixels. 

122 scienceBox : `lsst.images.Box` 

123 Bounding box of the science image. 

124 border : `float` 

125 Width, in input pixels, of the border of the input that warping left 

126 out of the coadd. 

127 

128 Returns 

129 ------- 

130 footprint : `lsst.images.Polygon` or `None` 

131 The part of ``scienceBox`` the input's warp covers, in science image 

132 pixels, or `None` if it covers none of it. 

133 """ 

134 inputBox = geom.Box2D(record.getBBox()) 

135 inputBox.grow(-border) 

136 inputBounds = Polygon.from_legacy(afwGeom.Polygon(inputBox)) 

137 polygon = record.getValidPolygon() 

138 inputToCoadd = _pixelTransform(record.getWcs(), coaddWcs) 

139 try: 

140 if polygon is not None: 140 ↛ 141line 140 didn't jump to line 141 because the condition on line 140 was never true

141 inputBounds = inputBounds.intersection(Polygon.from_legacy(polygon)) 

142 # Only vertices are transformed, so the edges of the input stay 

143 # straight. 

144 footprint = ( 

145 inputBounds.transform(inputToCoadd) 

146 .intersection(scienceFootprint) 

147 .transform(coaddToScienceTransform) 

148 ) 

149 # The convex hull keeps the footprint in one piece when it is clipped. 

150 footprint = footprint.convex_hull().intersection(scienceBox) 

151 except NoOverlapError: 

152 return None 

153 return Polygon.from_box(footprint) if isinstance(footprint, Box) else footprint 

154 

155 

156def _warpBorder(kernelName): 

157 """Return the width of the border of an image that warping leaves out. 

158 

159 A warped pixel is computed only where the whole kernel lies on the image. 

160 

161 Parameters 

162 ---------- 

163 kernelName : `str` 

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

165 

166 Returns 

167 ------- 

168 border : `float` 

169 Width of the border on each side, in pixels from the edge of the 

170 image's bounding box; 2.5 for ``lanczos3``. 

171 """ 

172 kernel = afwMath.WarpingControl(kernelName).getWarpingKernel() 

173 return kernel.getCtr().getX() + 0.5 

174 

175 

176@functools.cache 

177def _warpKernelSumSq(kernelName, nPhase=100): 

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

179 

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

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

182 diagonal became covariance between neighboring pixels, which a variance 

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

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

185 

186 Parameters 

187 ---------- 

188 kernelName : `str` 

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

190 nPhase : `int`, optional 

191 Number of sub-pixel phases to average over. 

192 

193 Returns 

194 ------- 

195 sumSq : `float` 

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

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

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

199 

200 Notes 

201 ----- 

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

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

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

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

206 """ 

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

208 # default cache size. 

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

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

211 total = 0.0 

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

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

214 kernel.computeImage(image, doNormalize=True) 

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

216 return (total/nPhase)**2 

217 

218 

219class GetTemplateConnections( 

220 pipeBase.PipelineTaskConnections, 

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

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

223): 

224 bbox = pipeBase.connectionTypes.Input( 

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

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

227 storageClass="Box2I", 

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

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

230 ) 

231 detector = pipeBase.connectionTypes.Input( 

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

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

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

235 storageClass="Detector", 

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

237 ) 

238 obs_info = pipeBase.connectionTypes.Input( 

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

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

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

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

243 storageClass="ObservationInfo", 

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

245 ) 

246 wcs = pipeBase.connectionTypes.Input( 

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

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

249 storageClass="Wcs", 

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

251 ) 

252 skyMap = pipeBase.connectionTypes.Input( 

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

254 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

255 dimensions=("skymap",), 

256 storageClass="SkyMap", 

257 ) 

258 coaddExposures = pipeBase.connectionTypes.Input( 

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

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

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

262 storageClass="ExposureF", 

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

264 multiple=True, 

265 deferLoad=True, 

266 deferGraphConstraint=True, 

267 ) 

268 

269 template = pipeBase.connectionTypes.Output( 

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

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

272 storageClass="ExposureF", 

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

274 ) 

275 

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

277 super().__init__(config=config) 

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

279 del self.bbox 

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

281 self.coaddExposures = dataclasses.replace( 

282 self.coaddExposures, 

283 deferGraphConstraint=False, 

284 ) 

285 if config.output_image_type == "future": 

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

287 else: 

288 del self.obs_info 

289 

290 

291class GetTemplateConfig( 

292 pipeBase.PipelineTaskConfig, pipelineConnections=GetTemplateConnections 

293): 

294 templateBorderSize = pexConfig.Field( 

295 dtype=int, 

296 default=20, 

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

298 ) 

299 warp = pexConfig.ConfigField( 

300 dtype=afwMath.Warper.ConfigClass, 

301 doc="warper configuration", 

302 ) 

303 doCorrectVarianceWarpCorrelation = pexConfig.Field( 

304 dtype=bool, 

305 default=True, 

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

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

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

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

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

311 " of the pixel scales involved.", 

312 ) 

313 coaddWarpKernel = pexConfig.ChoiceField( 

314 dtype=str, 

315 default="lanczos3", 

316 allowed={ 

317 "bilinear": "bilinear interpolation", 

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

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

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

321 }, 

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

323 " cannot be recovered from the coadds themselves. It is used to" 

324 " correct the variance if doCorrectVarianceWarpCorrelation is True," 

325 " and to trim the coadd inputs bboxes to account for the size of the" 

326 " warping kernel.", 

327 ) 

328 coaddPsf = pexConfig.ConfigField( 

329 doc="Configuration for CoaddPsf", 

330 dtype=CoaddPsfConfig, 

331 ) 

332 fwhmExposureGrid = pexConfig.Field( 

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

334 dtype=int, 

335 default=10, 

336 ) 

337 fwhmExposureBuffer = pexConfig.Field( 

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

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

340 dtype=float, 

341 default=0.05, 

342 ) 

343 varianceBackground = pexConfig.ConfigurableField( 

344 target=SubtractBackgroundTask, 

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

346 ) 

347 highVarianceThreshold = pexConfig.RangeField( 

348 dtype=float, 

349 default=4, 

350 min=1, 

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

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

353 ) 

354 highVarianceMaskFraction = pexConfig.Field( 

355 dtype=float, 

356 default=0.1, 

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

358 " HIGH_VARIANCE mask plane.", 

359 ) 

360 doScaleVariance = pexConfig.Field( 

361 dtype=bool, 

362 default=False, 

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

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

365 " doCorrectVarianceWarpCorrelation this overcounts the correction.", 

366 ) 

367 scaleVariance = pexConfig.ConfigurableField( 

368 target=ScaleVarianceTask, 

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

370 ) 

371 requireCoaddAtGraphBuild = pexConfig.Field( 

372 dtype=bool, 

373 default=False, 

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

375 " initial butler query during QuantumGraph generation.", 

376 ) 

377 output_image_type = pexConfig.ChoiceField[str]( 

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

379 allowed={ 

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

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

382 }, 

383 optional=False, 

384 default="legacy", 

385 ) 

386 raiseOnUndefinedMaskMap = pexConfig.Field( 

387 dtype=bool, 

388 default=False, 

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

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

391 " plane from the output instead." 

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

393 ) 

394 

395 def setDefaults(self): 

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

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

398 self.warp.cacheSize = 100000 

399 self.coaddPsf.cacheSize = self.warp.cacheSize 

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

401 # interpolation length for WCS evaluations. 

402 self.warp.interpLength = 100 

403 self.warp.warpingKernelName = "lanczos3" 

404 self.coaddPsf.warpingKernelName = self.warp.warpingKernelName 

405 

406 # Background subtraction of the variance plane 

407 self.varianceBackground.algorithm = "LINEAR" 

408 self.varianceBackground.binSize = 32 

409 self.varianceBackground.useApprox = False 

410 self.varianceBackground.statisticsProperty = "MEDIAN" 

411 self.varianceBackground.doFilterSuperPixels = True 

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

413 "EDGE", 

414 "DETECTED", 

415 "DETECTED_NEGATIVE", 

416 "NO_DATA", 

417 ] 

418 

419 

420class GetTemplateTask(pipeBase.PipelineTask): 

421 ConfigClass = GetTemplateConfig 

422 _DefaultName = "getTemplate" 

423 

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

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

426 if self.config.doScaleVariance: 

427 self.makeSubtask("scaleVariance") 

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

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

430 self.schema.addField( 

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

432 ) 

433 self.schema.addField( 

434 "patch", 

435 type=np.int32, 

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

437 ) 

438 self.schema.addField( 

439 "weight", 

440 type=float, 

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

442 ) 

443 self.schema.addField( 

444 "instrument", 

445 type=str, 

446 size=64, 

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

448 ) 

449 self.makeSubtask("varianceBackground") 

450 

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

452 inputs = butlerQC.get(inputRefs) 

453 detector = inputs.pop("detector") 

454 bbox = detector.getBBox() 

455 wcs = inputs.pop("wcs") 

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

457 coaddExposures = inputs.pop("coaddExposures") 

458 skymap = inputs.pop("skyMap") 

459 

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

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

462 

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

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

465 outputs = self.run( 

466 coaddExposureHandles=results.coaddExposures, 

467 bbox=bbox, 

468 wcs=wcs, 

469 dataIds=results.dataIds, 

470 physical_filter=physical_filter, 

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

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

473 ) 

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

475 dataId = butlerQC.quantum.dataId 

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

477 dataId.universe) 

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

479 inputRefs.coaddExposures, 

480 detector=detector, 

481 exposureRecord=exposureRecord, 

482 templateInputs=outputs.templateInputs) 

483 butlerQC.put(outputs, outputRefs) 

484 

485 def convert_outputs_to_future(self, template, coadd_refs, detector, exposureRecord, 

486 templateInputs=None): 

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

488 

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

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

491 instrument they were made from and their input image identifiers. 

492 The Gaussian PSF set by `run` is converted to an 

493 `lsst.images.psfs.GaussianPointSpreadFunction`. 

494 

495 Parameters 

496 ---------- 

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

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

499 and its mask is modified in place. 

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

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

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

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

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

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

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

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

508 for, which supplies the observation metadata recorded on the 

509 output. 

510 templateInputs : `dict` [`tuple` [`int`, `int`], \ 

511 `list` [`lsst.images.DifferenceImageTemplateInput`]], optional 

512 Input image records of each coadd, keyed by tract and patch, as 

513 returned by `run`. Coadds not listed record no inputs. 

514 

515 Returns 

516 ------- 

517 convertedTemplate : `lsst.images.DifferenceImage` 

518 The converted image. 

519 

520 Raises 

521 ------ 

522 RuntimeError 

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

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

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

526 spans more than one skymap. 

527 """ 

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

529 # template inputs to record which coadds contributed. 

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

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

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

533 instrument = catalog[0]["instrument"] 

534 template.setDetector(detector) 

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

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

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

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

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

540 for ref in coadd_refs} 

541 planeMap = get_legacy_template_mask_planes() 

542 if not self.config.raiseOnUndefinedMaskMap: 

543 self._clearUnmappedMaskPlanes(template.mask, planeMap) 

544 convertedTemplate = DifferenceImage.from_legacy( 

545 template, 

546 exposure_record=exposureRecord, 

547 unit=astropy.units.nJy, 

548 plane_map=planeMap, 

549 ) 

550 del template 

551 convertedTemplate.templates = DifferenceImageTemplateInfo.from_legacy_psf( 

552 convertedTemplate.sky_projection.pixel_frame, 

553 coaddPsf, 

554 skymaps.pop(), 

555 butler_info, 

556 log=self.log, 

557 instrument=instrument, 

558 ) 

559 if templateInputs: 

560 for info in convertedTemplate.templates: 

561 info.inputs = templateInputs.get((info.tract, info.patch), []) 

562 return convertedTemplate 

563 

564 def _clearUnmappedMaskPlanes(self, mask, planeMap): 

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

566 

567 Parameters 

568 ---------- 

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

570 Template mask to modify in place. 

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

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

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

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

575 """ 

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

577 for name in mask.getMaskPlaneDict(): 

578 if name in keep: 

579 continue 

580 bitmask = mask.getPlaneBitMask(name) 

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

582 if nSet: 

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

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

585 mask.array &= ~bitmask 

586 

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

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

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

590 structure of their dataIds. 

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

592 

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

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

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

596 all possibly-matching coadd exposures. 

597 

598 Parameters 

599 ---------- 

600 coaddExposureHandles : `iterable` \ 

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

602 `lsst.afw.image.Exposure`] 

603 Dataset handles to exposures that might overlap the desired 

604 region. 

605 bbox : `lsst.geom.Box2I` 

606 Template bounding box of the pixel geometry onto which the 

607 coaddExposures will be resampled. 

608 skymap : `lsst.skymap.SkyMap` 

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

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

611 Template WCS onto which the coadds will be resampled. 

612 

613 Returns 

614 ------- 

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

616 A struct with attributes: 

617 

618 ``coaddExposures`` 

619 Dict of coadd exposures that overlap the projected bbox, 

620 indexed on tract id 

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

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

623 ``dataIds`` 

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

625 projected bbox, indexed on tract id 

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

627 

628 Raises 

629 ------ 

630 NoWorkFound 

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

632 WCS is None. 

633 """ 

634 if wcs is None: 

635 raise pipeBase.NoWorkFound( 

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

637 ) 

638 

639 # Exposure's validPolygon would be more accurate 

640 detectorPolygon = geom.Box2D(bbox) 

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

642 overlappingArea = 0 

643 coaddExposures = collections.defaultdict(list) 

644 dataIds = collections.defaultdict(list) 

645 

646 for coaddRef in coaddExposureHandles: 

647 dataId = coaddRef.dataId 

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

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

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

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

652 # detector coordinates because the skymap's inverse mapping 

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

654 # arbitrary sky coordinates. 

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

656 if patchPolygon.intersection(detectorInPatchCoordinates): 

657 overlappingArea += patchPolygon.intersectionSingle( 

658 detectorInPatchCoordinates 

659 ).calculateArea() 

660 self.log.info( 

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

662 dataId["tract"], 

663 dataId["patch"], 

664 ) 

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

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

667 

668 if not overlappingArea: 

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

670 

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

672 

673 @timeMethod 

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

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

676 subtract from a science image. 

677 

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

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

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

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

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

683 CoaddPsf evaluated at the center of the detector. 

684 

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

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

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

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

689 unless none of the coadds records one. 

690 

691 Parameters 

692 ---------- 

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

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

695 `lsst.afw.image.Exposure`]] 

696 Coadds to be mosaicked, indexed on tract id. 

697 bbox : `lsst.geom.Box2I` 

698 Input Bounding box of the detector geometry onto which to 

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

700 will be grown by ``templateBorderSize``. 

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

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

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

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

705 tract id. 

706 physical_filter : `str` 

707 Physical filter of the science image. 

708 visit : `int`, optional 

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

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

711 science image for the visit. 

712 instrument : `str`, optional 

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

714 preference to coadds from any other instrument. 

715 

716 Returns 

717 ------- 

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

719 A struct with attributes: 

720 

721 ``template`` 

722 A template coadd exposure assembled out of patches, with 

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

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

725 ``templateInputs`` 

726 Input image records of each coadd that cover part of the 

727 science image, keyed by tract and patch (`dict` 

728 [`tuple` [`int`, `int`], 

729 `list` [`lsst.images.DifferenceImageTemplateInput`]]). Empty 

730 unless ``config.output_image_type`` is ``"future"``. 

731 

732 Raises 

733 ------ 

734 NoWorkFound 

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

736 RuntimeError 

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

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

739 images, so the template variance cannot be corrected. 

740 """ 

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

742 

743 scienceBBox = bbox 

744 bbox = geom.Box2I(bbox) 

745 bbox.grow(self.config.templateBorderSize) 

746 

747 warped = {} 

748 catalogs = [] 

749 inputs = [] 

750 plateScaleFactor = 1.0 if instrument == templateInstrument else None 

751 

752 for tract in coaddExposureHandles: 

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

754 coaddExposureHandles[tract], dataIds[tract], templateInstrument 

755 ) 

756 if len(catalog) == 0: 

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

758 continue 

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

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

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

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

763 if plateScaleFactor is None: 

764 plateScaleFactor = self._plateScaleFactor(coaddInputs[0], wcs, bbox) 

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

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

767 # Combine images from individual patches together. 

768 unwarped, count, included = self._merge( 

769 maskedImages, warpedBox, catalog[0].wcs 

770 ) 

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

772 del maskedImages 

773 if count == 0: 

774 self.log.info( 

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

776 tract, 

777 ) 

778 continue 

779 warpedBox.clip(totalBox) 

780 potentialInput = self.warper.warpExposure( 

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

782 ) 

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

784 del unwarped 

785 if np.all( 

786 potentialInput.mask.array 

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

788 ): 

789 self.log.info( 

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

791 tract, 

792 ) 

793 continue 

794 

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

796 tempCatalog = afwTable.ExposureCatalog(self.schema) 

797 tempCatalog.reserve(len(included)) 

798 for i in included: 

799 tempCatalog.append(catalog[i]) 

800 inputs.append(coaddInputs[i]) 

801 catalogs.append(tempCatalog) 

802 warped[tract] = potentialInput.maskedImage 

803 

804 if len(warped) == 0: 

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

806 

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

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

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

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

811 

812 if self.config.doScaleVariance: 

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

814 # needed. Note that the science variance is scaled 

815 # independently in ``AlardLuptonSubtractTask``. 

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

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

818 self.metadata["scaleTemplateVarianceFactor"] = varianceFactor 

819 

820 self._correctVariance(template, plateScaleFactor) 

821 

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

823 catalog = afwTable.ExposureCatalog(self.schema) 

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

825 for c in catalogs: 

826 catalog.extend(c) 

827 # Only the future output type records the inputs. 

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

829 templateInputs = self._collectTemplateInputs(catalog, inputs, wcs, scienceBBox) 

830 else: 

831 templateInputs = {} 

832 

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

834 self.checkHighVariance(template) 

835 if visit is not None: 

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

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

838 identity = afwImage.PhotoCalib(1.0) 

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

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

841 photoCalib = identity 

842 template.setPhotoCalib(photoCalib) 

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

844 template.setPsf(self._makeGaussianPsf(template)) 

845 

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

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

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

849 template.getInfo().setCoaddInputs(coaddInputs) 

850 return pipeBase.Struct(template=template, templateInputs=templateInputs) 

851 

852 def _correctVariance(self, template, plateScaleFactor): 

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

854 

855 Warping introduces a covariance between neighboring pixels, while the 

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

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

858 the detection significance on the convolved difference image is correct 

859 on average. 

860 

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

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

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

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

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

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

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

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

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

870 

871 Parameters 

872 ---------- 

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

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

875 plateScaleFactor : `float` or `None` 

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

877 

878 Raises 

879 ------ 

880 RuntimeError 

881 Raised if ``plateScaleFactor`` is `None`. 

882 """ 

883 if plateScaleFactor is None: 

884 raise RuntimeError( 

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

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

887 ) 

888 self.metadata["variancePlateScaleFactor"] = plateScaleFactor 

889 correlationFactor = 1.0 

890 if self.config.doCorrectVarianceWarpCorrelation: 

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

892 # and the one this task just performed. 

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

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

895 self.metadata["varianceWarpCorrelationFactor"] = correlationFactor 

896 scale = correlationFactor*plateScaleFactor 

897 if scale == 1.0: 

898 return 

899 template.variance.array *= scale 

900 self.metadata["templateVarianceCorrectionFactor"] = scale 

901 self.log.info( 

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

903 plateScaleFactor, 

904 ) 

905 self.log.info( 

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

907 correlationFactor, 

908 ) 

909 self.log.info( 

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

911 scale, 

912 ) 

913 

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

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

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

917 

918 Parameters 

919 ---------- 

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

921 Input records of one of the coadd patches. 

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

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