Coverage for python/lsst/pipe/tasks/computeExposureSummaryStats.py: 13%

1# This file is part of pipe_tasks. 

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/>. 

21 

22__all__ = ["ComputeExposureSummaryStatsTask", "ComputeExposureSummaryStatsConfig"] 

23 

24import warnings 

25import numpy as np 

26from scipy.stats import median_abs_deviation as sigmaMad 

27import astropy.units as units 

28from astropy.time import Time 

29from astropy.coordinates import AltAz, SkyCoord, EarthLocation 

30from lsst.daf.base import DateTime 

31 

32import lsst.pipe.base as pipeBase 

33import lsst.pex.config as pexConfig 

34import lsst.afw.math as afwMath 

35import lsst.afw.image as afwImage 

36import lsst.geom as geom 

37from lsst.meas.algorithms import ScienceSourceSelectorTask 

38from lsst.utils.timer import timeMethod 

39import lsst.ip.isr as ipIsr 

40 

41 

42class ComputeExposureSummaryStatsConfig(pexConfig.Config): 

43 """Config for ComputeExposureSummaryTask""" 

44 sigmaClip = pexConfig.Field( 

45 dtype=float, 

46 doc="Sigma for outlier rejection for sky noise.", 

47 default=3.0, 

48 ) 

49 clipIter = pexConfig.Field( 

50 dtype=int, 

51 doc="Number of iterations of outlier rejection for sky noise.", 

52 default=2, 

53 ) 

54 badMaskPlanes = pexConfig.ListField( 

55 dtype=str, 

56 doc="Mask planes that, if set, the associated pixel should not be included sky noise calculation.", 

57 default=("NO_DATA", "SUSPECT"), 

58 ) 

59 starSelection = pexConfig.Field( 

60 doc="Field to select full list of sources used for PSF modeling.", 

61 dtype=str, 

62 default="calib_psf_used", 

63 ) 

64 starSelector = pexConfig.ConfigurableField( 

65 target=ScienceSourceSelectorTask, 

66 doc="Selection of sources to compute PSF star statistics.", 

67 ) 

68 starShape = pexConfig.Field( 

69 doc="Base name of columns to use for the source shape in the PSF statistics computation.", 

70 dtype=str, 

71 default="slot_Shape" 

72 ) 

73 psfShape = pexConfig.Field( 

74 doc="Base name of columns to use for the PSF shape in the PSF statistics computation.", 

75 dtype=str, 

76 default="slot_PsfShape" 

77 ) 

78 psfSampling = pexConfig.Field( 

79 dtype=int, 

80 doc="Sampling rate in pixels in each dimension for the maxDistToNearestPsf metric " 

81 "caclulation grid (the tradeoff is between adequate sampling versus speed).", 

82 default=8, 

83 ) 

84 psfGridSampling = pexConfig.Field( 

85 dtype=int, 

86 doc="Sampling rate in pixels in each dimension for PSF model robustness metric " 

87 "caclulations grid (the tradeoff is between adequate sampling versus speed).", 

88 default=96, 

89 ) 

90 minPsfApRadiusPix = pexConfig.Field( 

91 dtype=float, 

92 doc="Minimum radius in pixels of the aperture within which to measure the flux of " 

93 "the PSF model for the psfApFluxDelta metric calculation (the radius is computed as " 

94 "max(``minPsfApRadius``, 3*psfSigma)).", 

95 default=2.0, 

96 ) 

97 psfApCorrFieldName = pexConfig.Field( 

98 doc="Name of the flux column associated with the aperture correction of the PSF model " 

99 "to use for the psfApCorrSigmaScaledDelta metric calculation.", 

100 dtype=str, 

101 default="base_PsfFlux_instFlux" 

102 ) 

103 psfBadMaskPlanes = pexConfig.ListField( 

104 dtype=str, 

105 doc="Mask planes that, if set, the associated pixel should not be included in the PSF model " 

106 "robutsness metric calculations (namely, maxDistToNearestPsf and psfTraceRadiusDelta).", 

107 default=("BAD", "CR", "EDGE", "INTRP", "NO_DATA", "SAT", "SUSPECT"), 

108 ) 

109 fiducialSkyBackground = pexConfig.DictField( 

110 keytype=str, 

111 itemtype=float, 

112 doc="Fiducial sky background level (ADU/s) assumed when calculating effective exposure time. " 

113 "Keyed by band.", 

114 default={'u': 1.0, 'g': 1.0, 'r': 1.0, 'i': 1.0, 'z': 1.0, 'y': 1.0}, 

115 ) 

116 fiducialPsfSigma = pexConfig.DictField( 

117 keytype=str, 

118 itemtype=float, 

119 doc="Fiducial PSF sigma (pixels) assumed when calculating effective exposure time. " 

120 "Keyed by band.", 

121 default={'u': 1.0, 'g': 1.0, 'r': 1.0, 'i': 1.0, 'z': 1.0, 'y': 1.0}, 

122 ) 

123 fiducialZeroPoint = pexConfig.DictField( 

124 keytype=str, 

125 itemtype=float, 

126 doc="Fiducial zero point assumed when calculating effective exposure time. " 

127 "Keyed by band.", 

128 default={'u': 25.0, 'g': 25.0, 'r': 25.0, 'i': 25.0, 'z': 25.0, 'y': 25.0}, 

129 ) 

130 maxEffectiveTransparency = pexConfig.Field( 

131 dtype=float, 

132 doc="Maximum value allowed for effective transparency scale factor (often inf or 1.0).", 

133 default=float('inf') 

134 ) 

135 magLimSnr = pexConfig.Field( 

136 dtype=float, 

137 doc="Signal-to-noise ratio for computing the magnitude limit depth.", 

138 default=5.0 

139 ) 

140 

141 def setDefaults(self): 

142 super().setDefaults() 

143 

144 self.starSelector.setDefaults() 

145 self.starSelector.doFlags = True 

146 self.starSelector.doSignalToNoise = True 

147 self.starSelector.doUnresolved = False 

148 self.starSelector.doIsolated = False 

149 self.starSelector.doRequireFiniteRaDec = False 

150 self.starSelector.doRequirePrimary = False 

151 

152 self.starSelector.signalToNoise.minimum = 50.0 

153 self.starSelector.signalToNoise.maximum = 1000.0 

154 

155 self.starSelector.flags.bad = ["slot_Shape_flag", "slot_PsfFlux_flag"] 

156 # Select stars used for PSF modeling. 

157 self.starSelector.flags.good = ["calib_psf_used"] 

158 

159 self.starSelector.signalToNoise.fluxField = "slot_PsfFlux_instFlux" 

160 self.starSelector.signalToNoise.errField = "slot_PsfFlux_instFluxErr" 

161 

162 

163class ComputeExposureSummaryStatsTask(pipeBase.Task): 

164 """Task to compute exposure summary statistics. 

165 

166 This task computes various quantities suitable for DPDD and other 

167 downstream processing at the detector centers, including: 

168 - expTime 

169 - psfSigma 

170 - psfArea 

171 - psfIxx 

172 - psfIyy 

173 - psfIxy 

174 - ra 

175 - dec 

176 - pixelScale (arcsec/pixel) 

177 - zenithDistance 

178 - zeroPoint 

179 - skyBg 

180 - skyNoise 

181 - meanVar 

182 - raCorners 

183 - decCorners 

184 - astromOffsetMean 

185 - astromOffsetStd 

186 

187 These additional quantities are computed from the stars in the detector: 

188 - psfStarDeltaE1Median 

189 - psfStarDeltaE2Median 

190 - psfStarDeltaE1Scatter 

191 - psfStarDeltaE2Scatter 

192 - psfStarDeltaSizeMedian 

193 - psfStarDeltaSizeScatter 

194 - psfStarScaledDeltaSizeScatter 

195 

196 These quantities are computed based on the PSF model and image mask 

197 to assess the robustness of the PSF model across a given detector 

198 (against, e.g., extrapolation instability): 

199 - maxDistToNearestPsf 

200 - psfTraceRadiusDelta 

201 - psfApFluxDelta 

202 

203 This quantity is computed based on the aperture correction map, the 

204 psfSigma, and the image mask to assess the robustness of the aperture 

205 corrections across a given detector: 

206 - psfApCorrSigmaScaledDelta 

207 

208 These quantities are computed to assess depth: 

209 - effTime 

210 - effTimePsfSigmaScale 

211 - effTimeSkyBgScale 

212 - effTimeZeroPointScale 

213 - magLim 

214 """ 

215 ConfigClass = ComputeExposureSummaryStatsConfig 

216 _DefaultName = "computeExposureSummaryStats" 

217 

218 def __init__(self, **kwargs): 

219 super().__init__(**kwargs) 

220 

221 self.makeSubtask("starSelector") 

222 

223 @timeMethod 

224 def run(self, exposure, sources, background): 

225 """Measure exposure statistics from the exposure, sources, and 

226 background. 

227 

228 Parameters 

229 ---------- 

230 exposure : `lsst.afw.image.ExposureF` 

231 sources : `lsst.afw.table.SourceCatalog` 

232 background : `lsst.afw.math.BackgroundList` 

233 

234 Returns 

235 ------- 

236 summary : `lsst.afw.image.ExposureSummary` 

237 """ 

238 self.log.info("Measuring exposure statistics") 

239 

240 summary = afwImage.ExposureSummaryStats() 

241 

242 # Set exposure time. 

243 exposureTime = exposure.getInfo().getVisitInfo().getExposureTime() 

244 summary.expTime = exposureTime 

245 

246 bbox = exposure.getBBox() 

247 

248 psf = exposure.getPsf() 

249 self.update_psf_stats( 

250 summary, psf, bbox, sources, image_mask=exposure.mask, image_ap_corr_map=exposure.apCorrMap 

251 ) 

252 

253 wcs = exposure.getWcs() 

254 visitInfo = exposure.getInfo().getVisitInfo() 

255 self.update_wcs_stats(summary, wcs, bbox, visitInfo) 

256 

257 photoCalib = exposure.getPhotoCalib() 

258 self.update_photo_calib_stats(summary, photoCalib) 

259 

260 self.update_background_stats(summary, background) 

261 

262 self.update_masked_image_stats(summary, exposure.getMaskedImage()) 

263 

264 self.update_magnitude_limit_stats(summary, exposure) 

265 

266 self.update_effective_time_stats(summary, exposure) 

267 

268 md = exposure.getMetadata() 

269 if 'SFM_ASTROM_OFFSET_MEAN' in md: 

270 summary.astromOffsetMean = md['SFM_ASTROM_OFFSET_MEAN'] 

271 summary.astromOffsetStd = md['SFM_ASTROM_OFFSET_STD'] 

272 

273 return summary 

274 

275 def update_psf_stats( 

276 self, 

277 summary, 

278 psf, 

279 bbox, 

280 sources=None, 

281 image_mask=None, 

282 image_ap_corr_map=None, 

283 sources_is_astropy=False, 

284 ): 

285 """Compute all summary-statistic fields that depend on the PSF model. 

286 

287 Parameters 

288 ---------- 

289 summary : `lsst.afw.image.ExposureSummaryStats` 

290 Summary object to update in-place. 

291 psf : `lsst.afw.detection.Psf` or `None` 

292 Point spread function model. If `None`, all fields that depend on 

293 the PSF will be reset (generally to NaN). 

294 bbox : `lsst.geom.Box2I` 

295 Bounding box of the image for which summary stats are being 

296 computed. 

297 sources : `lsst.afw.table.SourceCatalog` or `astropy.table.Table` 

298 Catalog for quantities that are computed from source table columns. 

299 If `None`, these quantities will be reset (generally to NaN). 

300 The type of this table must correspond to the 

301 ``sources_is_astropy`` argument. 

302 image_mask : `lsst.afw.image.Mask`, optional 

303 Mask image that may be used to compute distance-to-nearest-star 

304 metrics. 

305 sources_is_astropy : `bool`, optional 

306 Whether ``sources`` is an `astropy.table.Table` instance instead 

307 of an `lsst.afw.table.Catalog` instance. Default is `False` (the 

308 latter). 

309 """ 

310 nan = float("nan") 

311 summary.psfSigma = nan 

312 summary.psfIxx = nan 

313 summary.psfIyy = nan 

314 summary.psfIxy = nan 

315 summary.psfArea = nan 

316 summary.nPsfStar = 0 

317 summary.psfStarDeltaE1Median = nan 

318 summary.psfStarDeltaE2Median = nan 

319 summary.psfStarDeltaE1Scatter = nan 

320 summary.psfStarDeltaE2Scatter = nan 

321 summary.psfStarDeltaSizeMedian = nan 

322 summary.psfStarDeltaSizeScatter = nan 

323 summary.psfStarScaledDeltaSizeScatter = nan 

324 summary.maxDistToNearestPsf = nan 

325 summary.psfTraceRadiusDelta = nan 

326 summary.psfApFluxDelta = nan 

327 summary.psfApCorrSigmaScaledDelta = nan 

328 

329 if psf is None: 

330 return 

331 shape = psf.computeShape(bbox.getCenter()) 

332 summary.psfSigma = shape.getDeterminantRadius() 

333 summary.psfIxx = shape.getIxx() 

334 summary.psfIyy = shape.getIyy() 

335 summary.psfIxy = shape.getIxy() 

336 im = psf.computeKernelImage(bbox.getCenter()) 

337 # The calculation of effective psf area is taken from 

338 # ls.st/srd Equation 1. 

339 summary.psfArea = float(np.sum(im.array)**2./np.sum(im.array**2.)) 

340 

341 if image_mask is not None: 

342 psfApRadius = max(self.config.minPsfApRadiusPix, 3.0*summary.psfSigma) 

343 self.log.debug("Using radius of %.3f (pixels) for psfApFluxDelta metric", psfApRadius) 

344 psfTraceRadiusDelta, psfApFluxDelta = compute_psf_image_deltas( 

345 image_mask, 

346 psf, 

347 sampling=self.config.psfGridSampling, 

348 ap_radius_pix=psfApRadius, 

349 bad_mask_bits=self.config.psfBadMaskPlanes 

350 ) 

351 summary.psfTraceRadiusDelta = float(psfTraceRadiusDelta) 

352 summary.psfApFluxDelta = float(psfApFluxDelta) 

353 if image_ap_corr_map is not None: 

354 if self.config.psfApCorrFieldName not in image_ap_corr_map.keys(): 

355 self.log.warning(f"{self.config.psfApCorrFieldName} not found in " 

356 "image_ap_corr_map. Setting psfApCorrSigmaScaledDelta to NaN.") 

357 psfApCorrSigmaScaledDelta = nan 

358 else: 

359 image_ap_corr_field = image_ap_corr_map[self.config.psfApCorrFieldName] 

360 psfApCorrSigmaScaledDelta = compute_ap_corr_sigma_scaled_delta( 

361 image_mask, 

362 image_ap_corr_field, 

363 summary.psfSigma, 

364 sampling=self.config.psfGridSampling, 

365 bad_mask_bits=self.config.psfBadMaskPlanes, 

366 ) 

367 summary.psfApCorrSigmaScaledDelta = float(psfApCorrSigmaScaledDelta) 

368 

369 if sources is None: 

370 # No sources are available (as in some tests and rare cases where 

371 # the selection criteria in finalizeCharacterization lead to no 

372 # good sources). 

373 return 

374 

375 # Count the total number of psf stars used (prior to stats selection). 

376 nPsfStar = sources[self.config.starSelection].sum() 

377 summary.nPsfStar = int(nPsfStar) 

378 

379 psf_mask = self.starSelector.run(sources).selected 

380 nPsfStarsUsedInStats = psf_mask.sum() 

381 

382 if nPsfStarsUsedInStats == 0: 

383 # No stars to measure statistics, so we must return the defaults 

384 # of 0 stars and NaN values. 

385 return 

386 

387 if sources_is_astropy: 

388 psf_cat = sources[psf_mask] 

389 else: 

390 psf_cat = sources[psf_mask].copy(deep=True) 

391 

392 starXX = psf_cat[self.config.starShape + '_xx'] 

393 starYY = psf_cat[self.config.starShape + '_yy'] 

394 starXY = psf_cat[self.config.starShape + '_xy'] 

395 psfXX = psf_cat[self.config.psfShape + '_xx'] 

396 psfYY = psf_cat[self.config.psfShape + '_yy'] 

397 psfXY = psf_cat[self.config.psfShape + '_xy'] 

398 

399 # Use the trace radius for the star size. 

400 starSize = np.sqrt(starXX/2. + starYY/2.) 

401 

402 starE1 = (starXX - starYY)/(starXX + starYY) 

403 starE2 = 2*starXY/(starXX + starYY) 

404 starSizeMedian = np.median(starSize) 

405 

406 # Use the trace radius for the psf size. 

407 psfSize = np.sqrt(psfXX/2. + psfYY/2.) 

408 psfE1 = (psfXX - psfYY)/(psfXX + psfYY) 

409 psfE2 = 2*psfXY/(psfXX + psfYY) 

410 

411 psfStarDeltaE1Median = np.median(starE1 - psfE1) 

412 psfStarDeltaE1Scatter = sigmaMad(starE1 - psfE1, scale='normal') 

413 psfStarDeltaE2Median = np.median(starE2 - psfE2) 

414 psfStarDeltaE2Scatter = sigmaMad(starE2 - psfE2, scale='normal') 

415 

416 psfStarDeltaSizeMedian = np.median(starSize - psfSize) 

417 psfStarDeltaSizeScatter = sigmaMad(starSize - psfSize, scale='normal') 

418 psfStarScaledDeltaSizeScatter = psfStarDeltaSizeScatter/starSizeMedian 

419 

420 summary.psfStarDeltaE1Median = float(psfStarDeltaE1Median) 

421 summary.psfStarDeltaE2Median = float(psfStarDeltaE2Median) 

422 summary.psfStarDeltaE1Scatter = float(psfStarDeltaE1Scatter) 

423 summary.psfStarDeltaE2Scatter = float(psfStarDeltaE2Scatter) 

424 summary.psfStarDeltaSizeMedian = float(psfStarDeltaSizeMedian) 

425 summary.psfStarDeltaSizeScatter = float(psfStarDeltaSizeScatter) 

426 summary.psfStarScaledDeltaSizeScatter = float(psfStarScaledDeltaSizeScatter) 

427 

428 if image_mask is not None: 

429 maxDistToNearestPsf = maximum_nearest_psf_distance( 

430 image_mask, 

431 psf_cat, 

432 sampling=self.config.psfSampling, 

433 bad_mask_bits=self.config.psfBadMaskPlanes 

434 ) 

435 summary.maxDistToNearestPsf = float(maxDistToNearestPsf) 

436 

437 def update_wcs_stats(self, summary, wcs, bbox, visitInfo): 

438 """Compute all summary-statistic fields that depend on the WCS model. 

439 

440 Parameters 

441 ---------- 

442 summary : `lsst.afw.image.ExposureSummaryStats` 

443 Summary object to update in-place. 

444 wcs : `lsst.afw.geom.SkyWcs` or `None` 

445 Astrometric calibration model. If `None`, all fields that depend 

446 on the WCS will be reset (generally to NaN). 

447 bbox : `lsst.geom.Box2I` 

448 Bounding box of the image for which summary stats are being 

449 computed. 

450 visitInfo : `lsst.afw.image.VisitInfo` 

451 Observation information used in together with ``wcs`` to compute 

452 the zenith distance. 

453 """ 

454 nan = float("nan") 

455 summary.raCorners = [nan]*4 

456 summary.decCorners = [nan]*4 

457 summary.ra = nan 

458 summary.dec = nan 

459 summary.pixelScale = nan 

460 summary.zenithDistance = nan 

461 

462 if wcs is None: 

463 return 

464 

465 sph_pts = wcs.pixelToSky(geom.Box2D(bbox).getCorners()) 

466 summary.raCorners = [float(sph.getRa().asDegrees()) for sph in sph_pts] 

467 summary.decCorners = [float(sph.getDec().asDegrees()) for sph in sph_pts] 

468 

469 sph_pt = wcs.pixelToSky(bbox.getCenter()) 

470 summary.ra = sph_pt.getRa().asDegrees() 

471 summary.dec = sph_pt.getDec().asDegrees() 

472 summary.pixelScale = wcs.getPixelScale(bbox.getCenter()).asArcseconds() 

473 

474 date = visitInfo.getDate() 

475 

476 if date.isValid(): 

477 # We compute the zenithDistance at the center of the detector 

478 # rather than use the boresight value available via the visitInfo, 

479 # because the zenithDistance may vary significantly over a large 

480 # field of view. 

481 observatory = visitInfo.getObservatory() 

482 loc = EarthLocation(lat=observatory.getLatitude().asDegrees()*units.deg, 

483 lon=observatory.getLongitude().asDegrees()*units.deg, 

484 height=observatory.getElevation()*units.m) 

485 obstime = Time(visitInfo.getDate().get(system=DateTime.MJD), 

486 location=loc, format='mjd') 

487 coord = SkyCoord( 

488 summary.ra*units.degree, 

489 summary.dec*units.degree, 

490 obstime=obstime, 

491 location=loc, 

492 ) 

493 with warnings.catch_warnings(): 

494 warnings.simplefilter('ignore') 

495 altaz = coord.transform_to(AltAz) 

496 

497 summary.zenithDistance = float(90.0 - altaz.alt.degree) 

498 

499 def update_photo_calib_stats(self, summary, photo_calib): 

500 """Compute all summary-statistic fields that depend on the photometric 

501 calibration model. 

502 

503 Parameters 

504 ---------- 

505 summary : `lsst.afw.image.ExposureSummaryStats` 

506 Summary object to update in-place. 

507 photo_calib : `lsst.afw.image.PhotoCalib` or `None` 

508 Photometric calibration model. If `None`, all fields that depend 

509 on the photometric calibration will be reset (generally to NaN). 

510 """ 

511 if photo_calib is not None: 

512 summary.zeroPoint = float(2.5*np.log10(photo_calib.getInstFluxAtZeroMagnitude())) 

513 else: 

514 summary.zeroPoint = float("nan") 

515 

516 def update_background_stats(self, summary, background): 

517 """Compute summary-statistic fields that depend only on the 

518 background model. 

519 

520 Parameters 

521 ---------- 

522 summary : `lsst.afw.image.ExposureSummaryStats` 

523 Summary object to update in-place. 

524 background : `lsst.afw.math.BackgroundList` or `None` 

525 Background model. If `None`, all fields that depend on the 

526 background will be reset (generally to NaN). 

527 

528 Notes 

529 ----- 

530 This does not include fields that depend on the background-subtracted 

531 masked image; when the background changes, it should generally be 

532 applied to the image and `update_masked_image_stats` should be called 

533 as well. 

534 """ 

535 if background is not None: 

536 bgStats = (bg[0].getStatsImage().getImage().array 

537 for bg in background) 

538 summary.skyBg = float(sum(np.median(bg[np.isfinite(bg)]) for bg in bgStats)) 

539 else: 

540 summary.skyBg = float("nan") 

541 

542 def update_masked_image_stats(self, summary, masked_image): 

543 """Compute summary-statistic fields that depend on the masked image 

544 itself. 

545 

546 Parameters 

547 ---------- 

548 summary : `lsst.afw.image.ExposureSummaryStats` 

549 Summary object to update in-place. 

550 masked_image : `lsst.afw.image.MaskedImage` or `None` 

551 Masked image. If `None`, all fields that depend 

552 on the masked image will be reset (generally to NaN). 

553 """ 

554 nan = float("nan") 

555 if masked_image is None: 

556 summary.skyNoise = nan 

557 summary.meanVar = nan 

558 return 

559 statsCtrl = afwMath.StatisticsControl() 

560 statsCtrl.setNumSigmaClip(self.config.sigmaClip) 

561 statsCtrl.setNumIter(self.config.clipIter) 

562 statsCtrl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.badMaskPlanes)) 

563 statsCtrl.setNanSafe(True) 

564 

565 statObj = afwMath.makeStatistics(masked_image, afwMath.STDEVCLIP, statsCtrl) 

566 skyNoise, _ = statObj.getResult(afwMath.STDEVCLIP) 

567 summary.skyNoise = skyNoise 

568 

569 statObj = afwMath.makeStatistics(masked_image.variance, masked_image.mask, afwMath.MEANCLIP, 

570 statsCtrl) 

571 meanVar, _ = statObj.getResult(afwMath.MEANCLIP) 

572 summary.meanVar = meanVar 

573 

574 def update_effective_time_stats(self, summary, exposure): 

575 """Compute effective exposure time statistics to estimate depth. 

576 

577 The effective exposure time is the equivalent shutter open 

578 time that would be needed under nominal conditions to give the 

579 same signal-to-noise for a point source as what is achieved by 

580 the observation of interest. This metric combines measurements 

581 of the point-spread function, the sky brightness, and the 

582 transparency. It assumes that the observation is 

583 sky-background dominated. 

584 

585 .. _teff_definitions: 

586 

587 The effective exposure time and its subcomponents are defined in [1]_. 

588 

589 References 

590 ---------- 

591 

592 .. [1] Neilsen, E.H., Bernstein, G., Gruendl, R., and Kent, S. (2016). 

593 Limiting Magnitude, \tau, teff, and Image Quality in DES Year 1 

594 https://www.osti.gov/biblio/1250877/ 

595 

596 Parameters 

597 ---------- 

598 summary : `lsst.afw.image.ExposureSummaryStats` 

599 Summary object to update in-place. 

600 exposure : `lsst.afw.image.ExposureF` 

601 Exposure to grab band and exposure time metadata 

602 

603 """ 

604 nan = float("nan") 

605 summary.effTime = nan 

606 summary.effTimePsfSigmaScale = nan 

607 summary.effTimeSkyBgScale = nan 

608 summary.effTimeZeroPointScale = nan 

609 

610 exposureTime = exposure.getInfo().getVisitInfo().getExposureTime() 

611 filterLabel = exposure.getFilter() 

612 if (filterLabel is None) or (not filterLabel.hasBandLabel): 

613 band = None 

614 else: 

615 band = filterLabel.bandLabel 

616 

617 if band is None: 

618 self.log.warning("No band associated with exposure; effTime not calculated.") 

619 return 

620 

621 # PSF component 

622 if np.isnan(summary.psfSigma): 

623 self.log.debug("PSF sigma is NaN") 

624 f_eff = nan 

625 elif band not in self.config.fiducialPsfSigma: 

626 self.log.debug(f"Fiducial PSF value not found for {band}") 

627 f_eff = nan 

628 else: 

629 fiducialPsfSigma = self.config.fiducialPsfSigma[band] 

630 f_eff = (summary.psfSigma / fiducialPsfSigma)**-2 

631 

632 # Transparency component 

633 # Note: Assumes that the zeropoint includes the exposure time 

634 if np.isnan(summary.zeroPoint): 

635 self.log.debug("Zero point is NaN") 

636 c_eff = nan 

637 elif band not in self.config.fiducialZeroPoint: 

638 self.log.debug(f"Fiducial zero point value not found for {band}") 

639 c_eff = nan 

640 else: 

641 fiducialZeroPoint = self.config.fiducialZeroPoint[band] 

642 zeroPointDiff = fiducialZeroPoint - (summary.zeroPoint - 2.5*np.log10(exposureTime)) 

643 c_eff = min(10**(-2.0*(zeroPointDiff)/2.5), self.config.maxEffectiveTransparency) 

644 

645 # Sky brightness component (convert to cts/s) 

646 if np.isnan(summary.skyBg): 

647 self.log.debug("Sky background is NaN") 

648 b_eff = nan 

649 elif band not in self.config.fiducialSkyBackground: 

650 self.log.debug(f"Fiducial sky background value not found for {band}") 

651 b_eff = nan 

652 else: 

653 fiducialSkyBackground = self.config.fiducialSkyBackground[band] 

654 b_eff = fiducialSkyBackground/(summary.skyBg/exposureTime) 

655 

656 # Effective exposure time scale factor 

657 t_eff = f_eff * c_eff * b_eff 

658 

659 # Effective exposure time (seconds) 

660 effectiveTime = t_eff * exposureTime 

661 

662 # Output quantities 

663 summary.effTime = float(effectiveTime) 

664 summary.effTimePsfSigmaScale = float(f_eff) 

665 summary.effTimeSkyBgScale = float(b_eff) 

666 summary.effTimeZeroPointScale = float(c_eff) 

667 

668 def update_magnitude_limit_stats(self, summary, exposure): 

669 """Compute a summary statistic for the point-source magnitude limit at 

670 a given signal-to-noise ratio using exposure-level metadata. 

671 

672 The magnitude limit is calculated at a given SNR from the PSF, 

673 sky, zeropoint, and readnoise in accordance with SMTN-002 [1]_, 

674 LSE-40 [2]_, and DMTN-296 [3]_. 

675 

676 References 

677 ---------- 

678 

679 .. [1] "Calculating LSST limiting magnitudes and SNR" (SMTN-002) 

680 .. [2] "Photon Rates and SNR Calculations" (LSE-40) 

681 .. [3] "Calculations of Image and Catalog Depth" (DMTN-296) 

682 

683 

684 Parameters 

685 ---------- 

686 summary : `lsst.afw.image.ExposureSummaryStats` 

687 Summary object to update in-place. 

688 exposure : `lsst.afw.image.ExposureF` 

689 Exposure to grab band and exposure time metadata 

690 """ 

691 if exposure.getDetector() is None: 

692 summary.magLim = float("nan") 

693 return 

694 

695 # Calculate the average readnoise [e-] 

696 readNoiseList = list(ipIsr.getExposureReadNoises(exposure).values()) 

697 readNoise = np.nanmean(readNoiseList) 

698 if np.isnan(readNoise): 

699 readNoise = 0.0 

700 self.log.warning("Read noise set to NaN! Setting readNoise to 0.0 to proceed.") 

701 

702 # Calculate the average gain [e-/ADU] 

703 gainList = list(ipIsr.getExposureGains(exposure).values()) 

704 gain = np.nanmean(gainList) 

705 if np.isnan(gain): 

706 self.log.warning("Gain set to NaN! Setting magLim to NaN.") 

707 summary.magLim = float("nan") 

708 return 

709 

710 # Get the image units (default to 'adu' if metadata key absent) 

711 md = exposure.getMetadata() 

712 if md.get("LSST ISR UNITS", "adu") == "electron": 

713 gain = 1.0 

714 

715 # Convert readNoise to image units 

716 readNoise /= gain 

717 

718 # Calculate the limiting magnitude. 

719 # Note 1: Assumes that the image and readnoise have the same units 

720 # Note 2: Assumes that the zeropoint includes the exposure time 

721 magLim = compute_magnitude_limit(summary.psfArea, summary.skyBg, 

722 summary.zeroPoint, readNoise, 

723 gain, self.config.magLimSnr) 

724 

725 summary.magLim = float(magLim) 

726 

727 

728def maximum_nearest_psf_distance( 

729 image_mask, 

730 psf_cat, 

731 sampling=8, 

732 bad_mask_bits=["BAD", "CR", "INTRP", "SAT", "SUSPECT", "NO_DATA", "EDGE"], 

733): 

734 """Compute the maximum distance of an unmasked pixel to its nearest PSF. 

735 

736 Parameters 

737 ---------- 

738 image_mask : `lsst.afw.image.Mask` 

739 The mask plane associated with the exposure. 

740 psf_cat : `lsst.afw.table.SourceCatalog` or `astropy.table.Table` 

741 Catalog containing only the stars used in the PSF modeling. 

742 sampling : `int` 

743 Sampling rate in each dimension to create the grid of points on which 

744 to evaluate the distance to the nearest PSF star. The tradeoff is 

745 between adequate sampling versus speed. 

746 bad_mask_bits : `list` [`str`] 

747 Mask bits required to be absent for a pixel to be considered 

748 "unmasked". 

749 

750 Returns 

751 ------- 

752 max_dist_to_nearest_psf : `float` 

753 The maximum distance (in pixels) of an unmasked pixel to its nearest 

754 PSF model star. 

755 """ 

756 mask_arr = image_mask.array[::sampling, ::sampling] 

757 bitmask = image_mask.getPlaneBitMask(bad_mask_bits) 

758 good = ((mask_arr & bitmask) == 0) 

759 

760 x = np.arange(good.shape[1]) * sampling 

761 y = np.arange(good.shape[0]) * sampling 

762 xx, yy = np.meshgrid(x, y) 

763 

764 dist_to_nearest_psf = np.full(good.shape, np.inf) 

765 for psf in psf_cat: 

766 x_psf = psf["slot_Centroid_x"] 

767 y_psf = psf["slot_Centroid_y"] 

768 dist_to_nearest_psf = np.minimum(dist_to_nearest_psf, np.hypot(xx - x_psf, yy - y_psf)) 

769 unmasked_dists = dist_to_nearest_psf * good 

770 max_dist_to_nearest_psf = np.max(unmasked_dists) 

771 

772 return max_dist_to_nearest_psf 

773 

774 

775def compute_psf_image_deltas( 

776 image_mask, 

777 image_psf, 

778 sampling=96, 

779 ap_radius_pix=3.0, 

780 bad_mask_bits=["BAD", "CR", "INTRP", "SAT", "SUSPECT", "NO_DATA", "EDGE"], 

781): 

782 """Compute the delta between the maximum and minimum model PSF trace radius 

783 values evaluated on a grid of points lying in the unmasked region of the 

784 image. 

785 

786 Parameters 

787 ---------- 

788 image_mask : `lsst.afw.image.Mask` 

789 The mask plane associated with the exposure. 

790 image_psf : `lsst.afw.detection.Psf` 

791 The PSF model associated with the exposure. 

792 sampling : `int`, optional 

793 Sampling rate in each dimension to create the grid of points at which 

794 to evaluate ``image_psf``s trace radius value. The tradeoff is between 

795 adequate sampling versus speed. 

796 ap_radius_pix : `float`, optional 

797 Radius in pixels of the aperture on which to measure the flux of the 

798 PSF model. 

799 bad_mask_bits : `list` [`str`], optional 

800 Mask bits required to be absent for a pixel to be considered 

801 "unmasked". 

802 

803 Returns 

804 ------- 

805 psf_trace_radius_delta, psf_ap_flux_delta : `float` 

806 The delta (in pixels) between the maximum and minimum model PSF trace 

807 radius values and the PSF aperture fluxes (with aperture radius of 

808 max(2, 3*psfSigma)) evaluated on the x,y-grid subsampled on the 

809 unmasked detector pixels by a factor of ``sampling``. If both the 

810 model PSF trace radius value and aperture flux value on the grid 

811 evaluate to NaN, then NaNs are returned immediately. 

812 """ 

813 psf_trace_radius_list = [] 

814 psf_ap_flux_list = [] 

815 mask_arr = image_mask.array[::sampling, ::sampling] 

816 bitmask = image_mask.getPlaneBitMask(bad_mask_bits) 

817 good = ((mask_arr & bitmask) == 0) 

818 

819 x = np.arange(good.shape[1]) * sampling 

820 y = np.arange(good.shape[0]) * sampling 

821 xx, yy = np.meshgrid(x, y) 

822 

823 for x_mesh, y_mesh, good_mesh in zip(xx, yy, good): 

824 for x_point, y_point, is_good in zip(x_mesh, y_mesh, good_mesh): 

825 if is_good: 

826 position = geom.Point2D(x_point, y_point) 

827 psf_trace_radius = image_psf.computeShape(position).getTraceRadius() 

828 psf_ap_flux = image_psf.computeApertureFlux(ap_radius_pix, position) 

829 if ~np.isfinite(psf_trace_radius) and ~np.isfinite(psf_ap_flux): 

830 return float("nan"), float("nan") 

831 psf_trace_radius_list.append(psf_trace_radius) 

832 psf_ap_flux_list.append(psf_ap_flux) 

833 

834 psf_trace_radius_delta = np.max(psf_trace_radius_list) - np.min(psf_trace_radius_list) 

835 if np.any(np.asarray(psf_ap_flux_list) < 0.0): # Consider any -ve flux entry as "bad". 

836 psf_ap_flux_delta = float("nan") 

837 else: 

838 psf_ap_flux_delta = np.max(psf_ap_flux_list) - np.min(psf_ap_flux_list) 

839 

840 return psf_trace_radius_delta, psf_ap_flux_delta 

841 

842 

843def compute_ap_corr_sigma_scaled_delta( 

844 image_mask, 

845 image_ap_corr_field, 

846 psfSigma, 

847 sampling=96, 

848 bad_mask_bits=["BAD", "CR", "INTRP", "SAT", "SUSPECT", "NO_DATA", "EDGE"], 

849): 

850 """Compute the delta between the maximum and minimum aperture correction 

851 values scaled (divided) by ``psfSigma`` for the given field representation, 

852 ``image_ap_corr_field`` evaluated on a grid of points lying in the 

853 unmasked region of the image. 

854 

855 Parameters 

856 ---------- 

857 image_mask : `lsst.afw.image.Mask` 

858 The mask plane associated with the exposure. 

859 image_ap_corr_field : `lsst.afw.math.ChebyshevBoundedField` 

860 The ChebyshevBoundedField representation of the aperture correction 

861 of interest for the exposure. 

862 psfSigma : `float` 

863 The PSF model second-moments determinant radius (center of chip) 

864 in pixels. 

865 sampling : `int`, optional 

866 Sampling rate in each dimension to create the grid of points at which 

867 to evaluate ``image_psf``s trace radius value. The tradeoff is between 

868 adequate sampling versus speed. 

869 bad_mask_bits : `list` [`str`], optional 

870 Mask bits required to be absent for a pixel to be considered 

871 "unmasked". 

872 

873 Returns