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 fiducialReadNoise = pexConfig.DictField( 

131 keytype=str, 

132 itemtype=float, 

133 doc="Fiducial readnoise (electrons) assumed when calculating effective exposure time. " 

134 "Keyed by band.", 

135 default={'u': 9.0, 'g': 9.0, 'r': 9.0, 'i': 9.0, 'z': 9.0, 'y': 9.0}, 

136 ) 

137 fiducialExpTime = pexConfig.DictField( 

138 keytype=str, 

139 itemtype=float, 

140 doc="Fiducial exposure time (seconds). " 

141 "Keyed by band.", 

142 default={'u': 30.0, 'g': 30.0, 'r': 30.0, 'i': 30.0, 'z': 30.0, 'y': 30.0}, 

143 ) 

144 fiducialMagLim = pexConfig.DictField( 

145 keytype=str, 

146 itemtype=float, 

147 doc="Fiducial magnitude limit depth at SNR=5. " 

148 "Keyed by band.", 

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

150 ) 

151 maxEffectiveTransparency = pexConfig.Field( 

152 dtype=float, 

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

154 default=float('inf') 

155 ) 

156 magLimSnr = pexConfig.Field( 

157 dtype=float, 

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

159 default=5.0 

160 ) 

161 

162 def setDefaults(self): 

163 super().setDefaults() 

164 

165 self.starSelector.setDefaults() 

166 self.starSelector.doFlags = True 

167 self.starSelector.doSignalToNoise = True 

168 self.starSelector.doUnresolved = False 

169 self.starSelector.doIsolated = False 

170 self.starSelector.doRequireFiniteRaDec = False 

171 self.starSelector.doRequirePrimary = False 

172 

173 self.starSelector.signalToNoise.minimum = 50.0 

174 self.starSelector.signalToNoise.maximum = 1000.0 

175 

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

177 # Select stars used for PSF modeling. 

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

179 

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

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

182 

183 def fiducialMagnitudeLimit(self, band, pixelScale, gain): 

184 """Compute the fiducial point-source magnitude limit based on config values. 

185 This follows the conventions laid out in SMTN-002, LSE-40, and DMTN-296. 

186 

187 Parameters 

188 ---------- 

189 band : `str` 

190 The band of interest 

191 pixelScale : `float` 

192 The pixel scale [arcsec/pix] 

193 gain : `float` 

194 The instrumental gain for the exposure [e-/ADU]. The gain should 

195 be 1.0 if the image units are [e-]. 

196 

197 Returns 

198 ------- 

199 magnitude_limit : `float` 

200 The fiducial magnitude limit calculated from fiducial values. 

201 """ 

202 nan = float("nan") 

203 

204 # Fiducials from config 

205 fiducialPsfSigma = self.fiducialPsfSigma.get(band, nan) 

206 fiducialSkyBackground = self.fiducialSkyBackground.get(band, nan) 

207 fiducialZeroPoint = self.fiducialZeroPoint.get(band, nan) 

208 fiducialReadNoise = self.fiducialReadNoise.get(band, nan) 

209 fiducialExpTime = self.fiducialExpTime.get(band, nan) 

210 magLimSnr = self.magLimSnr 

211 

212 # Derived fiducial quantities 

213 fiducialPsfArea = psf_sigma_to_psf_area(fiducialPsfSigma, pixelScale) 

214 fiducialZeroPoint += 2.5*np.log10(fiducialExpTime) 

215 fiducialSkyBg = fiducialSkyBackground * fiducialExpTime 

216 fiducialReadNoise /= gain 

217 

218 # Calculate the fiducial magnitude limit 

219 fiducialMagLim = compute_magnitude_limit(fiducialPsfArea, 

220 fiducialSkyBg, 

221 fiducialZeroPoint, 

222 fiducialReadNoise, 

223 gain, 

224 magLimSnr) 

225 

226 return fiducialMagLim 

227 

228 

229class ComputeExposureSummaryStatsTask(pipeBase.Task): 

230 """Task to compute exposure summary statistics. 

231 

232 This task computes various quantities suitable for DPDD and other 

233 downstream processing at the detector centers, including: 

234 - expTime 

235 - psfSigma 

236 - psfArea 

237 - psfIxx 

238 - psfIyy 

239 - psfIxy 

240 - ra 

241 - dec 

242 - pixelScale (arcsec/pixel) 

243 - zenithDistance 

244 - zeroPoint 

245 - skyBg 

246 - skyNoise 

247 - meanVar 

248 - raCorners 

249 - decCorners 

250 - astromOffsetMean 

251 - astromOffsetStd 

252 

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

254 - psfStarDeltaE1Median 

255 - psfStarDeltaE2Median 

256 - psfStarDeltaE1Scatter 

257 - psfStarDeltaE2Scatter 

258 - psfStarDeltaSizeMedian 

259 - psfStarDeltaSizeScatter 

260 - psfStarScaledDeltaSizeScatter 

261 

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

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

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

265 - maxDistToNearestPsf 

266 - psfTraceRadiusDelta 

267 - psfApFluxDelta 

268 

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

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

271 corrections across a given detector: 

272 - psfApCorrSigmaScaledDelta 

273 

274 These quantities are computed based on the shape measurements of the 

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

276 measure of the deviation from a circular shape): 

277 - starEMedian 

278 - starUnNormalizedEMedian 

279 

280 These quantities are computed to assess depth: 

281 - effTime 

282 - effTimePsfSigmaScale 

283 - effTimeSkyBgScale 

284 - effTimeZeroPointScale 

285 - magLim 

286 """ 

287 ConfigClass = ComputeExposureSummaryStatsConfig 

288 _DefaultName = "computeExposureSummaryStats" 

289 

290 def __init__(self, **kwargs): 

291 super().__init__(**kwargs) 

292 

293 self.makeSubtask("starSelector") 

294 

295 @timeMethod 

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

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

298 background. 

299 

300 Parameters 

301 ---------- 

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

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

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

305 

306 Returns 

307 ------- 

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

309 """ 

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

311 

312 summary = afwImage.ExposureSummaryStats() 

313 

314 # Set exposure time. 

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

316 summary.expTime = exposureTime 

317 

318 bbox = exposure.getBBox() 

319 

320 psf = exposure.getPsf() 

321 self.update_psf_stats( 

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

323 ) 

324 

325 wcs = exposure.getWcs() 

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

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

328 

329 photoCalib = exposure.getPhotoCalib() 

330 self.update_photo_calib_stats(summary, photoCalib) 

331 

332 self.update_background_stats(summary, background) 

333 

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

335 

336 self.update_magnitude_limit_stats(summary, exposure) 

337 

338 self.update_effective_time_stats(summary, exposure) 

339 

340 md = exposure.getMetadata() 

341 if 'SFM_ASTROM_OFFSET_MEAN' in md: 

342 summary.astromOffsetMean = md['SFM_ASTROM_OFFSET_MEAN'] 

343 summary.astromOffsetStd = md['SFM_ASTROM_OFFSET_STD'] 

344 

345 return summary 

346 

347 def update_psf_stats( 

348 self, 

349 summary, 

350 psf, 

351 bbox, 

352 sources=None, 

353 image_mask=None, 

354 image_ap_corr_map=None, 

355 sources_is_astropy=False, 

356 ): 

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

358 

359 Parameters 

360 ---------- 

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

362 Summary object to update in-place. 

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

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

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

366 bbox : `lsst.geom.Box2I` 

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

368 computed. 

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

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

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

372 The type of this table must correspond to the 

373 ``sources_is_astropy`` argument. 

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

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

376 metrics. 

377 sources_is_astropy : `bool`, optional 

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

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

380 latter). 

381 """ 

382 nan = float("nan") 

383 summary.psfSigma = nan 

384 summary.psfIxx = nan 

385 summary.psfIyy = nan 

386 summary.psfIxy = nan 

387 summary.psfArea = nan 

388 summary.nPsfStar = 0 

389 summary.psfStarDeltaE1Median = nan 

390 summary.psfStarDeltaE2Median = nan 

391 summary.psfStarDeltaE1Scatter = nan 

392 summary.psfStarDeltaE2Scatter = nan 

393 summary.psfStarDeltaSizeMedian = nan 

394 summary.psfStarDeltaSizeScatter = nan 

395 summary.psfStarScaledDeltaSizeScatter = nan 

396 summary.maxDistToNearestPsf = nan 

397 summary.psfTraceRadiusDelta = nan 

398 summary.psfApFluxDelta = nan 

399 summary.psfApCorrSigmaScaledDelta = nan 

400 summary.starEMedian = nan 

401 summary.starUnNormalizedEMedian = nan 

402 

403 if psf is None: 

404 return 

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

406 summary.psfSigma = shape.getDeterminantRadius() 

407 summary.psfIxx = shape.getIxx() 

408 summary.psfIyy = shape.getIyy() 

409 summary.psfIxy = shape.getIxy() 

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

411 # The calculation of effective psf area is taken from 

412 # ls.st/srd Equation 1. 

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

414 

415 if image_mask is not None: 

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

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

418 psfTraceRadiusDelta, psfApFluxDelta = compute_psf_image_deltas( 

419 image_mask, 

420 psf, 

421 sampling=self.config.psfGridSampling, 

422 ap_radius_pix=psfApRadius, 

423 bad_mask_bits=self.config.psfBadMaskPlanes 

424 ) 

425 summary.psfTraceRadiusDelta = float(psfTraceRadiusDelta) 

426 summary.psfApFluxDelta = float(psfApFluxDelta) 

427 if image_ap_corr_map is not None: 

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

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

430 "image_ap_corr_map. Setting psfApCorrSigmaScaledDelta to NaN.") 

431 psfApCorrSigmaScaledDelta = nan 

432 else: 

433 image_ap_corr_field = image_ap_corr_map[self.config.psfApCorrFieldName] 

434 psfApCorrSigmaScaledDelta = compute_ap_corr_sigma_scaled_delta( 

435 image_mask, 

436 image_ap_corr_field, 

437 summary.psfSigma, 

438 sampling=self.config.psfGridSampling, 

439 bad_mask_bits=self.config.psfBadMaskPlanes, 

440 ) 

441 summary.psfApCorrSigmaScaledDelta = float(psfApCorrSigmaScaledDelta) 

442 

443 if sources is None: 

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

445 # the selection criteria in finalizeCharacterization lead to no 

446 # good sources). 

447 return 

448 

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

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

451 summary.nPsfStar = int(nPsfStar) 

452 

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

454 nPsfStarsUsedInStats = psf_mask.sum() 

455 

456 if nPsfStarsUsedInStats == 0: 

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

458 # of 0 stars and NaN values. 

459 return 

460 

461 if sources_is_astropy: 

462 psf_cat = sources[psf_mask] 

463 else: 

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

465 

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

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

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

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

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

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

472 

473 # Use the trace radius for the star size. 

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

475 

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

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

478 starSizeMedian = np.median(starSize) 

479 

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

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

482 starEMedian = np.median(starE) 

483 starUnNormalizedEMedian = np.median(starUnNormalizedE) 

484 summary.starEMedian = float(starEMedian) 

485 summary.starUnNormalizedEMedian = float(starUnNormalizedEMedian) 

486 

487 # Use the trace radius for the psf size. 

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

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

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

491 

492 psfStarDeltaE1Median = np.median(starE1 - psfE1) 

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

494 psfStarDeltaE2Median = np.median(starE2 - psfE2) 

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

496 

497 psfStarDeltaSizeMedian = np.median(starSize - psfSize) 

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

499 psfStarScaledDeltaSizeScatter = psfStarDeltaSizeScatter/starSizeMedian 

500 

501 summary.psfStarDeltaE1Median = float(psfStarDeltaE1Median) 

502 summary.psfStarDeltaE2Median = float(psfStarDeltaE2Median) 

503 summary.psfStarDeltaE1Scatter = float(psfStarDeltaE1Scatter) 

504 summary.psfStarDeltaE2Scatter = float(psfStarDeltaE2Scatter) 

505 summary.psfStarDeltaSizeMedian = float(psfStarDeltaSizeMedian) 

506 summary.psfStarDeltaSizeScatter = float(psfStarDeltaSizeScatter) 

507 summary.psfStarScaledDeltaSizeScatter = float(psfStarScaledDeltaSizeScatter) 

508 

509 if image_mask is not None: 

510 maxDistToNearestPsf = maximum_nearest_psf_distance( 

511 image_mask, 

512 psf_cat, 

513 sampling=self.config.psfSampling, 

514 bad_mask_bits=self.config.psfBadMaskPlanes 

515 ) 

516 summary.maxDistToNearestPsf = float(maxDistToNearestPsf) 

517 

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

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

520 

521 Parameters 

522 ---------- 

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

524 Summary object to update in-place. 

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

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

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

528 bbox : `lsst.geom.Box2I` 

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

530 computed. 

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

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

533 the zenith distance. 

534 """ 

535 nan = float("nan") 

536 summary.raCorners = [nan]*4 

537 summary.decCorners = [nan]*4 

538 summary.ra = nan 

539 summary.dec = nan 

540 summary.pixelScale = nan 

541 summary.zenithDistance = nan 

542 

543 if wcs is None: 

544 return 

545 

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

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

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

549 

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

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

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

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

554 

555 date = visitInfo.getDate() 

556 

557 if date.isValid(): 

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

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

560 # because the zenithDistance may vary significantly over a large 

561 # field of view. 

562 observatory = visitInfo.getObservatory() 

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

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

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

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

567 location=loc, format='mjd') 

568 coord = SkyCoord( 

569 summary.ra*units.degree, 

570 summary.dec*units.degree, 

571 obstime=obstime, 

572 location=loc, 

573 ) 

574 with warnings.catch_warnings(): 

575 warnings.simplefilter('ignore') 

576 altaz = coord.transform_to(AltAz) 

577 

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

579 

580 def update_photo_calib_stats(self, summary, photo_calib): 

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

582 calibration model. 

583 

584 Parameters 

585 ---------- 

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

587 Summary object to update in-place. 

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

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

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

591 """ 

592 if photo_calib is not None: 

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

594 else: 

595 summary.zeroPoint = float("nan") 

596 

597 def update_background_stats(self, summary, background): 

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

599 background model. 

600 

601 Parameters 

602 ---------- 

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

604 Summary object to update in-place. 

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

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

607 background will be reset (generally to NaN). 

608 

609 Notes 

610 ----- 

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

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

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

614 as well. 

615 """ 

616 if background is not None: 

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

618 for bg in background) 

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

620 else: 

621 summary.skyBg = float("nan") 

622 

623 def update_masked_image_stats(self, summary, masked_image): 

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

625 itself. 

626 

627 Parameters 

628 ---------- 

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

630 Summary object to update in-place. 

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

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

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

634 """ 

635 nan = float("nan") 

636 if masked_image is None: 

637 summary.skyNoise = nan 

638 summary.meanVar = nan 

639 return 

640 statsCtrl = afwMath.StatisticsControl() 

641 statsCtrl.setNumSigmaClip(self.config.sigmaClip) 

642 statsCtrl.setNumIter(self.config.clipIter) 

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

644 statsCtrl.setNanSafe(True) 

645 

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

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

648 summary.skyNoise = skyNoise 

649 

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

651 statsCtrl) 

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

653 summary.meanVar = meanVar 

654 

655 def update_effective_time_stats(self, summary, exposure): 

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

657 

658 The effective exposure time is the equivalent shutter open 

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

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

661 the observation of interest. This metric combines measurements 

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

663 transparency. It assumes that the observation is 

664 sky-background dominated. 

665 

666 .. _teff_definitions: 

667 

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

669 

670 References 

671 ---------- 

672 

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

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

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

676 

677 Parameters 

678 ---------- 

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

680 Summary object to update in-place. 

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

682 Exposure to grab band and exposure time metadata 

683 

684 """ 

685 nan = float("nan") 

686 summary.effTime = nan 

687 summary.effTimePsfSigmaScale = nan 

688 summary.effTimeSkyBgScale = nan 

689 summary.effTimeZeroPointScale = nan 

690 

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

692 filterLabel = exposure.getFilter() 

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

694 band = None 

695 else: 

696 band = filterLabel.bandLabel 

697 

698 if band is None: 

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

700 return 

701 

702 # PSF component 

703 if np.isnan(summary.psfSigma): 

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

705 f_eff = nan 

706 elif band not in self.config.fiducialPsfSigma: 

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

708 f_eff = nan 

709 else: 

710 fiducialPsfSigma = self.config.fiducialPsfSigma[band] 

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

712 

713 # Transparency component 

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

715 if np.isnan(summary.zeroPoint): 

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

717 c_eff = nan 

718 elif band not in self.config.fiducialZeroPoint: 

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

720 c_eff = nan 

721 else: 

722 fiducialZeroPoint = self.config.fiducialZeroPoint[band] 

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

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

725 

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

727 if np.isnan(summary.skyBg): 

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

729 b_eff = nan 

730 elif band not in self.config.fiducialSkyBackground: 

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

732 b_eff = nan 

733 else: 

734 fiducialSkyBackground = self.config.fiducialSkyBackground[band] 

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

736 

737 # Old effective exposure time scale factor 

738 # t_eff = f_eff * c_eff * b_eff 

739 

740 # New effective exposure time (seconds) from magnitude limit 

741 if band not in self.config.fiducialMagLim: 

742 self.log.debug(f"Fiducial magnitude limit not found for {band}") 

743 effectiveTime = nan 

744 elif band not in self.config.fiducialExpTime: 

745 self.log.debug(f"Fiducial exposure time not found for {band}") 

746 effectiveTime = nan 

747 else: 

748 effectiveTime = compute_effective_time(summary.magLim, 

749 self.config.fiducialMagLim[band], 

750 self.config.fiducialExpTime[band]) 

751 

752 # Output quantities 

753 summary.effTime = float(effectiveTime) 

754 summary.effTimePsfSigmaScale = float(f_eff) 

755 summary.effTimeSkyBgScale = float(b_eff) 

756 summary.effTimeZeroPointScale = float(c_eff) 

757 

758 def update_magnitude_limit_stats(self, summary, exposure): 

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

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

761 

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

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

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

765 

766 References 

767 ---------- 

768 

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

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

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

772 

773 

774 Parameters 

775 ---------- 

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

777 Summary object to update in-place. 

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

779 Exposure to grab band and exposure time metadata 

780 """ 

781 if exposure.getDetector() is None: 

782 summary.magLim = float("nan") 

783 return 

784 

785 # Calculate the average readnoise [e-] 

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

787 readNoise = np.nanmean(readNoiseList) 

788 if np.isnan(readNoise): 

789 readNoise = 0.0 

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

791 

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

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

794 gain = np.nanmean(gainList) 

795 if np.isnan(gain): 

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

797 summary.magLim = float("nan") 

798 return 

799 

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

801 md = exposure.getMetadata() 

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

803 gain = 1.0 

804 

805 # Convert readNoise to image units 

806 readNoise /= gain 

807 

808 # Calculate the limiting magnitude. 

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

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

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

812 summary.zeroPoint, readNoise, 

813 gain, self.config.magLimSnr) 

814 

815 summary.magLim = float(magLim) 

816 

817 

818def maximum_nearest_psf_distance( 

819 image_mask, 

820 psf_cat, 

821 sampling=8, 

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

823): 

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

825 

826 Parameters 

827 ---------- 

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

829 The mask plane associated with the exposure. 

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

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

832 sampling : `int` 

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

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

835 between adequate sampling versus speed. 

836 bad_mask_bits : `list` [`str`] 

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

838 "unmasked". 

839 

840 Returns 

841 ------- 

842 max_dist_to_nearest_psf : `float` 

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

844 PSF model star. 

845 """ 

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

847 bitmask = image_mask.getPlaneBitMask(bad_mask_bits) 

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

849 

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

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

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

853 

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

855 for psf in psf_cat: 

856 x_psf = psf["slot_Centroid_x"] 

857 y_psf = psf["slot_Centroid_y"] 

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

859 unmasked_dists = dist_to_nearest_psf * good