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 based on the shape measurements of the 

209 sources used in the PSF model to assess the image quality (i.e. as a 

210 measure of the deviation from a circular shape): 

211 - starEMedian 

212 - starUnNormalizedEMedian 

213 

214 These quantities are computed to assess depth: 

215 - effTime 

216 - effTimePsfSigmaScale 

217 - effTimeSkyBgScale 

218 - effTimeZeroPointScale 

219 - magLim 

220 """ 

221 ConfigClass = ComputeExposureSummaryStatsConfig 

222 _DefaultName = "computeExposureSummaryStats" 

223 

224 def __init__(self, **kwargs): 

225 super().__init__(**kwargs) 

226 

227 self.makeSubtask("starSelector") 

228 

229 @timeMethod 

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

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

232 background. 

233 

234 Parameters 

235 ---------- 

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

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

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

239 

240 Returns 

241 ------- 

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

243 """ 

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

245 

246 summary = afwImage.ExposureSummaryStats() 

247 

248 # Set exposure time. 

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

250 summary.expTime = exposureTime 

251 

252 bbox = exposure.getBBox() 

253 

254 psf = exposure.getPsf() 

255 self.update_psf_stats( 

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

257 ) 

258 

259 wcs = exposure.getWcs() 

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

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

262 

263 photoCalib = exposure.getPhotoCalib() 

264 self.update_photo_calib_stats(summary, photoCalib) 

265 

266 self.update_background_stats(summary, background) 

267 

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

269 

270 self.update_magnitude_limit_stats(summary, exposure) 

271 

272 self.update_effective_time_stats(summary, exposure) 

273 

274 md = exposure.getMetadata() 

275 if 'SFM_ASTROM_OFFSET_MEAN' in md: 

276 summary.astromOffsetMean = md['SFM_ASTROM_OFFSET_MEAN'] 

277 summary.astromOffsetStd = md['SFM_ASTROM_OFFSET_STD'] 

278 

279 return summary 

280 

281 def update_psf_stats( 

282 self, 

283 summary, 

284 psf, 

285 bbox, 

286 sources=None, 

287 image_mask=None, 

288 image_ap_corr_map=None, 

289 sources_is_astropy=False, 

290 ): 

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

292 

293 Parameters 

294 ---------- 

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

296 Summary object to update in-place. 

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

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

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

300 bbox : `lsst.geom.Box2I` 

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

302 computed. 

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

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

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

306 The type of this table must correspond to the 

307 ``sources_is_astropy`` argument. 

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

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

310 metrics. 

311 sources_is_astropy : `bool`, optional 

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

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

314 latter). 

315 """ 

316 nan = float("nan") 

317 summary.psfSigma = nan 

318 summary.psfIxx = nan 

319 summary.psfIyy = nan 

320 summary.psfIxy = nan 

321 summary.psfArea = nan 

322 summary.nPsfStar = 0 

323 summary.psfStarDeltaE1Median = nan 

324 summary.psfStarDeltaE2Median = nan 

325 summary.psfStarDeltaE1Scatter = nan 

326 summary.psfStarDeltaE2Scatter = nan 

327 summary.psfStarDeltaSizeMedian = nan 

328 summary.psfStarDeltaSizeScatter = nan 

329 summary.psfStarScaledDeltaSizeScatter = nan 

330 summary.maxDistToNearestPsf = nan 

331 summary.psfTraceRadiusDelta = nan 

332 summary.psfApFluxDelta = nan 

333 summary.psfApCorrSigmaScaledDelta = nan 

334 summary.starEMedian = nan 

335 summary.starUnNormalizedEMedian = nan 

336 

337 if psf is None: 

338 return 

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

340 summary.psfSigma = shape.getDeterminantRadius() 

341 summary.psfIxx = shape.getIxx() 

342 summary.psfIyy = shape.getIyy() 

343 summary.psfIxy = shape.getIxy() 

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

345 # The calculation of effective psf area is taken from 

346 # ls.st/srd Equation 1. 

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

348 

349 if image_mask is not None: 

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

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

352 psfTraceRadiusDelta, psfApFluxDelta = compute_psf_image_deltas( 

353 image_mask, 

354 psf, 

355 sampling=self.config.psfGridSampling, 

356 ap_radius_pix=psfApRadius, 

357 bad_mask_bits=self.config.psfBadMaskPlanes 

358 ) 

359 summary.psfTraceRadiusDelta = float(psfTraceRadiusDelta) 

360 summary.psfApFluxDelta = float(psfApFluxDelta) 

361 if image_ap_corr_map is not None: 

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

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

364 "image_ap_corr_map. Setting psfApCorrSigmaScaledDelta to NaN.") 

365 psfApCorrSigmaScaledDelta = nan 

366 else: 

367 image_ap_corr_field = image_ap_corr_map[self.config.psfApCorrFieldName] 

368 psfApCorrSigmaScaledDelta = compute_ap_corr_sigma_scaled_delta( 

369 image_mask, 

370 image_ap_corr_field, 

371 summary.psfSigma, 

372 sampling=self.config.psfGridSampling, 

373 bad_mask_bits=self.config.psfBadMaskPlanes, 

374 ) 

375 summary.psfApCorrSigmaScaledDelta = float(psfApCorrSigmaScaledDelta) 

376 

377 if sources is None: 

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

379 # the selection criteria in finalizeCharacterization lead to no 

380 # good sources). 

381 return 

382 

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

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

385 summary.nPsfStar = int(nPsfStar) 

386 

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

388 nPsfStarsUsedInStats = psf_mask.sum() 

389 

390 if nPsfStarsUsedInStats == 0: 

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

392 # of 0 stars and NaN values. 

393 return 

394 

395 if sources_is_astropy: 

396 psf_cat = sources[psf_mask] 

397 else: 

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

399 

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

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

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

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

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

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

406 

407 # Use the trace radius for the star size. 

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

409 

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

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

412 starSizeMedian = np.median(starSize) 

413 

414 starE = np.sqrt(starE1**2.0 + starE2**2.0) 

415 starUnNormalizedE = np.sqrt((starXX - starYY)**2.0 + (2.0*starXY)**2.0) 

416 starEMedian = np.median(starE) 

417 starUnNormalizedEMedian = np.median(starUnNormalizedE) 

418 summary.starEMedian = float(starEMedian) 

419 summary.starUnNormalizedEMedian = float(starUnNormalizedEMedian) 

420 

421 # Use the trace radius for the psf size. 

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

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

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

425 

426 psfStarDeltaE1Median = np.median(starE1 - psfE1) 

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

428 psfStarDeltaE2Median = np.median(starE2 - psfE2) 

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

430 

431 psfStarDeltaSizeMedian = np.median(starSize - psfSize) 

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

433 psfStarScaledDeltaSizeScatter = psfStarDeltaSizeScatter/starSizeMedian 

434 

435 summary.psfStarDeltaE1Median = float(psfStarDeltaE1Median) 

436 summary.psfStarDeltaE2Median = float(psfStarDeltaE2Median) 

437 summary.psfStarDeltaE1Scatter = float(psfStarDeltaE1Scatter) 

438 summary.psfStarDeltaE2Scatter = float(psfStarDeltaE2Scatter) 

439 summary.psfStarDeltaSizeMedian = float(psfStarDeltaSizeMedian) 

440 summary.psfStarDeltaSizeScatter = float(psfStarDeltaSizeScatter) 

441 summary.psfStarScaledDeltaSizeScatter = float(psfStarScaledDeltaSizeScatter) 

442 

443 if image_mask is not None: 

444 maxDistToNearestPsf = maximum_nearest_psf_distance( 

445 image_mask, 

446 psf_cat, 

447 sampling=self.config.psfSampling, 

448 bad_mask_bits=self.config.psfBadMaskPlanes 

449 ) 

450 summary.maxDistToNearestPsf = float(maxDistToNearestPsf) 

451 

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

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

454 

455 Parameters 

456 ---------- 

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

458 Summary object to update in-place. 

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

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

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

462 bbox : `lsst.geom.Box2I` 

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

464 computed. 

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

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

467 the zenith distance. 

468 """ 

469 nan = float("nan") 

470 summary.raCorners = [nan]*4 

471 summary.decCorners = [nan]*4 

472 summary.ra = nan 

473 summary.dec = nan 

474 summary.pixelScale = nan 

475 summary.zenithDistance = nan 

476 

477 if wcs is None: 

478 return 

479 

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

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

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

483 

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

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

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

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

488 

489 date = visitInfo.getDate() 

490 

491 if date.isValid(): 

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

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

494 # because the zenithDistance may vary significantly over a large 

495 # field of view. 

496 observatory = visitInfo.getObservatory() 

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

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

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

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

501 location=loc, format='mjd') 

502 coord = SkyCoord( 

503 summary.ra*units.degree, 

504 summary.dec*units.degree, 

505 obstime=obstime, 

506 location=loc, 

507 ) 

508 with warnings.catch_warnings(): 

509 warnings.simplefilter('ignore') 

510 altaz = coord.transform_to(AltAz) 

511 

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

513 

514 def update_photo_calib_stats(self, summary, photo_calib): 

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

516 calibration model. 

517 

518 Parameters 

519 ---------- 

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

521 Summary object to update in-place. 

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

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

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

525 """ 

526 if photo_calib is not None: 

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

528 else: 

529 summary.zeroPoint = float("nan") 

530 

531 def update_background_stats(self, summary, background): 

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

533 background model. 

534 

535 Parameters 

536 ---------- 

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

538 Summary object to update in-place. 

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

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

541 background will be reset (generally to NaN). 

542 

543 Notes 

544 ----- 

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

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

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

548 as well. 

549 """ 

550 if background is not None: 

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

552 for bg in background) 

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

554 else: 

555 summary.skyBg = float("nan") 

556 

557 def update_masked_image_stats(self, summary, masked_image): 

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

559 itself. 

560 

561 Parameters 

562 ---------- 

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

564 Summary object to update in-place. 

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

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

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

568 """ 

569 nan = float("nan") 

570 if masked_image is None: 

571 summary.skyNoise = nan 

572 summary.meanVar = nan 

573 return 

574 statsCtrl = afwMath.StatisticsControl() 

575 statsCtrl.setNumSigmaClip(self.config.sigmaClip) 

576 statsCtrl.setNumIter(self.config.clipIter) 

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

578 statsCtrl.setNanSafe(True) 

579 

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

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

582 summary.skyNoise = skyNoise 

583 

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

585 statsCtrl) 

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

587 summary.meanVar = meanVar 

588 

589 def update_effective_time_stats(self, summary, exposure): 

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

591 

592 The effective exposure time is the equivalent shutter open 

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

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

595 the observation of interest. This metric combines measurements 

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

597 transparency. It assumes that the observation is 

598 sky-background dominated. 

599 

600 .. _teff_definitions: 

601 

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

603 

604 References 

605 ---------- 

606 

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

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

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

610 

611 Parameters 

612 ---------- 

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

614 Summary object to update in-place. 

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

616 Exposure to grab band and exposure time metadata 

617 

618 """ 

619 nan = float("nan") 

620 summary.effTime = nan 

621 summary.effTimePsfSigmaScale = nan 

622 summary.effTimeSkyBgScale = nan 

623 summary.effTimeZeroPointScale = nan 

624 

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

626 filterLabel = exposure.getFilter() 

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

628 band = None 

629 else: 

630 band = filterLabel.bandLabel 

631 

632 if band is None: 

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

634 return 

635 

636 # PSF component 

637 if np.isnan(summary.psfSigma): 

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

639 f_eff = nan 

640 elif band not in self.config.fiducialPsfSigma: 

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

642 f_eff = nan 

643 else: 

644 fiducialPsfSigma = self.config.fiducialPsfSigma[band] 

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

646 

647 # Transparency component 

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

649 if np.isnan(summary.zeroPoint): 

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

651 c_eff = nan 

652 elif band not in self.config.fiducialZeroPoint: 

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

654 c_eff = nan 

655 else: 

656 fiducialZeroPoint = self.config.fiducialZeroPoint[band] 

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

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

659 

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

661 if np.isnan(summary.skyBg): 

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

663 b_eff = nan 

664 elif band not in self.config.fiducialSkyBackground: 

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

666 b_eff = nan 

667 else: 

668 fiducialSkyBackground = self.config.fiducialSkyBackground[band] 

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

670 

671 # Effective exposure time scale factor 

672 t_eff = f_eff * c_eff * b_eff 

673 

674 # Effective exposure time (seconds) 

675 effectiveTime = t_eff * exposureTime 

676 

677 # Output quantities 

678 summary.effTime = float(effectiveTime) 

679 summary.effTimePsfSigmaScale = float(f_eff) 

680 summary.effTimeSkyBgScale = float(b_eff) 

681 summary.effTimeZeroPointScale = float(c_eff) 

682 

683 def update_magnitude_limit_stats(self, summary, exposure): 

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

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

686 

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

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

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

690 

691 References 

692 ---------- 

693 

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

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

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

697 

698 

699 Parameters 

700 ---------- 

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

702 Summary object to update in-place. 

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

704 Exposure to grab band and exposure time metadata 

705 """ 

706 if exposure.getDetector() is None: 

707 summary.magLim = float("nan") 

708 return 

709 

710 # Calculate the average readnoise [e-] 

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

712 readNoise = np.nanmean(readNoiseList) 

713 if np.isnan(readNoise): 

714 readNoise = 0.0 

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

716 

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

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

719 gain = np.nanmean(gainList) 

720 if np.isnan(gain): 

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

722 summary.magLim = float("nan") 

723 return 

724 

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

726 md = exposure.getMetadata() 

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

728 gain = 1.0 

729 

730 # Convert readNoise to image units 

731 readNoise /= gain 

732 

733 # Calculate the limiting magnitude. 

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

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

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

737 summary.zeroPoint, readNoise, 

738 gain, self.config.magLimSnr) 

739 

740 summary.magLim = float(magLim) 

741 

742 

743def maximum_nearest_psf_distance( 

744 image_mask, 

745 psf_cat, 

746 sampling=8, 

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

748): 

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

750 

751 Parameters 

752 ---------- 

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

754 The mask plane associated with the exposure. 

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

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

757 sampling : `int` 

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

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

760 between adequate sampling versus speed. 

761 bad_mask_bits : `list` [`str`] 

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

763 "unmasked". 

764 

765 Returns 

766 ------- 

767 max_dist_to_nearest_psf : `float` 

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

769 PSF model star. 

770 """ 

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

772 bitmask = image_mask.getPlaneBitMask(bad_mask_bits) 

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

774 

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

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

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

778 

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

780 for psf in psf_cat: 

781 x_psf = psf["slot_Centroid_x"] 

782 y_psf = psf["slot_Centroid_y"] 

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

784 unmasked_dists = dist_to_nearest_psf * good 

785 max_dist_to_nearest_psf = np.max(unmasked_dists) 

786 

787 return max_dist_to_nearest_psf 

788 

789 

790def compute_psf_image_deltas( 

791 image_mask, 

792 image_psf, 

793 sampling=96, 

794 ap_radius_pix=3.0, 

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

796): 

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

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

799 image. 

800 

801 Parameters 

802 ---------- 

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

804 The mask plane associated with the exposure. 

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

806 The PSF model associated with the exposure. 

807 sampling : `int`, optional 

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

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

810 adequate sampling versus speed. 

811 ap_radius_pix : `float`, optional 

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

813 PSF model. 

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

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

816 "unmasked". 

817 

818 Returns 

819 ------- 

820 psf_trace_radius_delta, psf_ap_flux_delta : `float` 

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

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

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

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

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

826 evaluate to NaN, then NaNs are returned immediately. 

827 """ 

828 psf_trace_radius_list = [] 

829 psf_ap_flux_list = [] 

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

831 bitmask = image_mask.getPlaneBitMask(bad_mask_bits) 

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

833 

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

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

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

837 

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

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

840 if is_good: 

841 position = geom.Point2D(x_point, y_point) 

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

843 psf_ap_flux = image_psf.computeApertureFlux(ap_radius_pix, position) 

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

845 return float("nan"), float("nan") 

846 psf_trace_radius_list.append(psf_trace_radius) 

847 psf_ap_flux_list.append(psf_ap_flux) 

848 

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

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

851 psf_ap_flux_delta = float("nan") 

852 else: 

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

854 

855 return psf_trace_radius_delta, psf_ap_flux_delta 

856 

857 

858def compute_ap_corr_sigma_scaled_delta( 

859 image_mask, 

860 image_ap_corr_field, 

861 psfSigma, 

862 sampling=96, 

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