Coverage for python / lsst / pipe / tasks / computeExposureSummaryStats.py: 12%
412 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-05 08:30 +0000
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/>.
22__all__ = ["ComputeExposureSummaryStatsTask", "ComputeExposureSummaryStatsConfig"]
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
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
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 starHigherOrderMomentBase = pexConfig.Field(
79 doc="Base name of the columns from which to obtain the higher-order moments.",
80 dtype=str,
81 default="ext_shapeHSM_HigherOrderMomentsSource",
82 )
83 psfSampling = pexConfig.Field(
84 dtype=int,
85 doc="Sampling rate in pixels in each dimension for the maxDistToNearestPsf metric "
86 "caclulation grid (the tradeoff is between adequate sampling versus speed).",
87 default=8,
88 )
89 psfGridSampling = pexConfig.Field(
90 dtype=int,
91 doc="Sampling rate in pixels in each dimension for PSF model robustness metric "
92 "caclulations grid (the tradeoff is between adequate sampling versus speed).",
93 default=96,
94 )
95 minPsfApRadiusPix = pexConfig.Field(
96 dtype=float,
97 doc="Minimum radius in pixels of the aperture within which to measure the flux of "
98 "the PSF model for the psfApFluxDelta metric calculation (the radius is computed as "
99 "max(``minPsfApRadius``, 3*psfSigma)).",
100 default=2.0,
101 )
102 psfApCorrFieldName = pexConfig.Field(
103 doc="Name of the flux column associated with the aperture correction of the PSF model "
104 "to use for the psfApCorrSigmaScaledDelta metric calculation.",
105 dtype=str,
106 default="base_PsfFlux_instFlux"
107 )
108 psfBadMaskPlanes = pexConfig.ListField(
109 dtype=str,
110 doc="Mask planes that, if set, the associated pixel should not be included in the PSF model "
111 "robutsness metric calculations (namely, maxDistToNearestPsf and psfTraceRadiusDelta).",
112 default=("BAD", "CR", "EDGE", "INTRP", "NO_DATA", "SAT", "SUSPECT"),
113 )
114 fiducialSkyBackground = pexConfig.DictField(
115 keytype=str,
116 itemtype=float,
117 doc="Fiducial sky background level (ADU/s) assumed when calculating effective exposure time. "
118 "Keyed by band.",
119 default={"u": 1.0, "g": 1.0, "r": 1.0, "i": 1.0, "z": 1.0, "y": 1.0},
120 )
121 fiducialPsfSigma = pexConfig.DictField(
122 keytype=str,
123 itemtype=float,
124 doc="Fiducial PSF sigma (pixels) assumed when calculating effective exposure time. "
125 "Keyed by band.",
126 default={"u": 1.0, "g": 1.0, "r": 1.0, "i": 1.0, "z": 1.0, "y": 1.0},
127 )
128 fiducialZeroPoint = pexConfig.DictField(
129 keytype=str,
130 itemtype=float,
131 doc="Fiducial zero point assumed when calculating effective exposure time. "
132 "Keyed by band.",
133 default={"u": 25.0, "g": 25.0, "r": 25.0, "i": 25.0, "z": 25.0, "y": 25.0},
134 )
135 fiducialReadNoise = pexConfig.DictField(
136 keytype=str,
137 itemtype=float,
138 doc="Fiducial readnoise (electrons) assumed when calculating effective exposure time. "
139 "Keyed by band.",
140 default={"u": 9.0, "g": 9.0, "r": 9.0, "i": 9.0, "z": 9.0, "y": 9.0},
141 )
142 fiducialExpTime = pexConfig.DictField(
143 keytype=str,
144 itemtype=float,
145 doc="Fiducial exposure time (seconds). "
146 "Keyed by band.",
147 default={"u": 30.0, "g": 30.0, "r": 30.0, "i": 30.0, "z": 30.0, "y": 30.0},
148 )
149 fiducialMagLim = pexConfig.DictField(
150 keytype=str,
151 itemtype=float,
152 doc="Fiducial magnitude limit depth at SNR=5. "
153 "Keyed by band.",
154 default={"u": 25.0, "g": 25.0, "r": 25.0, "i": 25.0, "z": 25.0, "y": 25.0},
155 )
156 maxEffectiveTransparency = pexConfig.Field(
157 dtype=float,
158 doc="Maximum value allowed for effective transparency scale factor (often inf or 1.0).",
159 default=float("inf")
160 )
161 magLimSnr = pexConfig.Field(
162 dtype=float,
163 doc="Signal-to-noise ratio for computing the magnitude limit depth.",
164 default=5.0
165 )
167 def setDefaults(self):
168 super().setDefaults()
170 self.starSelector.setDefaults()
171 self.starSelector.doFlags = True
172 self.starSelector.doSignalToNoise = True
173 self.starSelector.doUnresolved = False
174 self.starSelector.doIsolated = False
175 self.starSelector.doRequireFiniteRaDec = False
176 self.starSelector.doRequirePrimary = False
178 self.starSelector.signalToNoise.minimum = 50.0
179 self.starSelector.signalToNoise.maximum = 1000.0
181 self.starSelector.flags.bad = ["slot_Shape_flag", "slot_PsfFlux_flag"]
182 # Select stars used for PSF modeling.
183 self.starSelector.flags.good = ["calib_psf_used"]
185 self.starSelector.signalToNoise.fluxField = "slot_PsfFlux_instFlux"
186 self.starSelector.signalToNoise.errField = "slot_PsfFlux_instFluxErr"
188 def fiducialMagnitudeLimit(self, band, pixelScale, gain):
189 """Compute the fiducial point-source magnitude limit based on config values.
190 This follows the conventions laid out in SMTN-002, LSE-40, and DMTN-296.
192 Parameters
193 ----------
194 band : `str`
195 The band of interest
196 pixelScale : `float`
197 The pixel scale [arcsec/pix]
198 gain : `float`
199 The instrumental gain for the exposure [e-/ADU]. The gain should
200 be 1.0 if the image units are [e-].
202 Returns
203 -------
204 magnitude_limit : `float`
205 The fiducial magnitude limit calculated from fiducial values.
206 """
207 nan = float("nan")
209 # Fiducials from config
210 fiducialPsfSigma = self.fiducialPsfSigma.get(band, nan)
211 fiducialSkyBackground = self.fiducialSkyBackground.get(band, nan)
212 fiducialZeroPoint = self.fiducialZeroPoint.get(band, nan)
213 fiducialReadNoise = self.fiducialReadNoise.get(band, nan)
214 fiducialExpTime = self.fiducialExpTime.get(band, nan)
215 magLimSnr = self.magLimSnr
217 # Derived fiducial quantities
218 fiducialPsfArea = psf_sigma_to_psf_area(fiducialPsfSigma, pixelScale)
219 fiducialZeroPoint += 2.5*np.log10(fiducialExpTime)
220 fiducialSkyBg = fiducialSkyBackground * fiducialExpTime
221 fiducialReadNoise /= gain
223 # Calculate the fiducial magnitude limit
224 fiducialMagLim = compute_magnitude_limit(fiducialPsfArea,
225 fiducialSkyBg,
226 fiducialZeroPoint,
227 fiducialReadNoise,
228 gain,
229 magLimSnr)
231 return fiducialMagLim
234class ComputeExposureSummaryStatsTask(pipeBase.Task):
235 """Task to compute exposure summary statistics.
237 This task computes various quantities suitable for DPDD and other
238 downstream processing at the detector centers, including:
239 - expTime
240 - psfSigma
241 - psfArea
242 - psfIxx
243 - psfIyy
244 - psfIxy
245 - ra
246 - dec
247 - pixelScale (arcsec/pixel)
248 - zenithDistance
249 - zeroPoint
250 - skyBg
251 - skyNoise
252 - meanVar
253 - raCorners
254 - decCorners
255 - astromOffsetMean
256 - astromOffsetStd
258 These additional quantities are computed from the stars in the detector:
259 - psfStarDeltaE1Median
260 - psfStarDeltaE2Median
261 - psfStarDeltaE1Scatter
262 - psfStarDeltaE2Scatter
263 - psfStarDeltaSizeMedian
264 - psfStarDeltaSizeScatter
265 - psfStarScaledDeltaSizeScatter
267 These quantities are computed based on the PSF model and image mask
268 to assess the robustness of the PSF model across a given detector
269 (against, e.g., extrapolation instability):
270 - maxDistToNearestPsf
271 - psfTraceRadiusDelta
272 - psfApFluxDelta
274 This quantity is computed based on the aperture correction map, the
275 psfSigma, and the image mask to assess the robustness of the aperture
276 corrections across a given detector:
277 - psfApCorrSigmaScaledDelta
279 These quantities are computed based on the shape measurements of the
280 sources used in the PSF model to assess the image quality (i.e. as a
281 measure of the deviation from a circular shape):
282 - starEMedian
283 - starUnNormalizedEMedian
285 These quantities are computed to assess depth:
286 - effTime
287 - effTimePsfSigmaScale
288 - effTimeSkyBgScale
289 - effTimeZeroPointScale
290 - magLim
291 """
292 ConfigClass = ComputeExposureSummaryStatsConfig
293 _DefaultName = "computeExposureSummaryStats"
295 def __init__(self, **kwargs):
296 super().__init__(**kwargs)
298 self.makeSubtask("starSelector")
300 @timeMethod
301 def run(self, exposure, sources, background):
302 """Measure exposure statistics from the exposure, sources, and
303 background.
305 Parameters
306 ----------
307 exposure : `lsst.afw.image.ExposureF`
308 sources : `lsst.afw.table.SourceCatalog`
309 background : `lsst.afw.math.BackgroundList`
311 Returns
312 -------
313 summary : `lsst.afw.image.ExposureSummary`
314 """
315 self.log.info("Measuring exposure statistics")
317 summary = afwImage.ExposureSummaryStats()
319 # Set exposure time.
320 exposureTime = exposure.getInfo().getVisitInfo().getExposureTime()
321 summary.expTime = exposureTime
323 bbox = exposure.getBBox()
325 psf = exposure.getPsf()
326 self.update_psf_stats(
327 summary, psf, bbox, sources, image_mask=exposure.mask, image_ap_corr_map=exposure.apCorrMap
328 )
330 wcs = exposure.getWcs()
331 visitInfo = exposure.getInfo().getVisitInfo()
332 self.update_wcs_stats(summary, wcs, bbox, visitInfo)
334 photoCalib = exposure.getPhotoCalib()
335 self.update_photo_calib_stats(summary, photoCalib)
337 self.update_background_stats(summary, background)
339 self.update_masked_image_stats(summary, exposure.getMaskedImage())
341 self.update_magnitude_limit_stats(summary, exposure)
343 self.update_effective_time_stats(summary, exposure)
345 md = exposure.getMetadata()
346 if "SFM_ASTROM_OFFSET_MEAN" in md:
347 summary.astromOffsetMean = md["SFM_ASTROM_OFFSET_MEAN"]
348 summary.astromOffsetStd = md["SFM_ASTROM_OFFSET_STD"]
350 return summary
352 def update_psf_stats(
353 self,
354 summary,
355 psf,
356 bbox,
357 sources=None,
358 image_mask=None,
359 image_ap_corr_map=None,
360 sources_is_astropy=False,
361 ):
362 """Compute all summary-statistic fields that depend on the PSF model.
364 Parameters
365 ----------
366 summary : `lsst.afw.image.ExposureSummaryStats`
367 Summary object to update in-place.
368 psf : `lsst.afw.detection.Psf` or `None`
369 Point spread function model. If `None`, all fields that depend on
370 the PSF will be reset (generally to NaN).
371 bbox : `lsst.geom.Box2I`
372 Bounding box of the image for which summary stats are being
373 computed.
374 sources : `lsst.afw.table.SourceCatalog` or `astropy.table.Table`
375 Catalog for quantities that are computed from source table columns.
376 If `None`, these quantities will be reset (generally to NaN).
377 The type of this table must correspond to the
378 ``sources_is_astropy`` argument.
379 image_mask : `lsst.afw.image.Mask`, optional
380 Mask image that may be used to compute distance-to-nearest-star
381 metrics.
382 sources_is_astropy : `bool`, optional
383 Whether ``sources`` is an `astropy.table.Table` instance instead
384 of an `lsst.afw.table.Catalog` instance. Default is `False` (the
385 latter).
386 """
387 nan = float("nan")
388 summary.psfSigma = nan
389 summary.psfIxx = nan
390 summary.psfIyy = nan
391 summary.psfIxy = nan
392 summary.psfArea = nan
393 summary.nPsfStar = 0
394 summary.psfStarDeltaE1Median = nan
395 summary.psfStarDeltaE2Median = nan
396 summary.psfStarDeltaE1Scatter = nan
397 summary.psfStarDeltaE2Scatter = nan
398 summary.psfStarDeltaSizeMedian = nan
399 summary.psfStarDeltaSizeScatter = nan
400 summary.psfStarScaledDeltaSizeScatter = nan
401 summary.maxDistToNearestPsf = nan
402 summary.psfTraceRadiusDelta = nan
403 summary.psfApFluxDelta = nan
404 summary.psfApCorrSigmaScaledDelta = nan
405 summary.starEMedian = nan
406 summary.starUnNormalizedEMedian = nan
407 summary.starComa1Median = nan
408 summary.starComa2Median = nan
409 summary.starTrefoil1Median = nan
410 summary.starTrefoil2Median = nan
411 summary.starKurtosisMedian = nan
412 summary.starE41Median = nan
413 summary.starE42Median = nan
415 if psf is None:
416 return
417 shape = psf.computeShape(bbox.getCenter())
418 summary.psfSigma = shape.getDeterminantRadius()
419 summary.psfIxx = shape.getIxx()
420 summary.psfIyy = shape.getIyy()
421 summary.psfIxy = shape.getIxy()
422 im = psf.computeKernelImage(bbox.getCenter())
423 # The calculation of effective psf area is taken from
424 # ls.st/srd Equation 1.
425 summary.psfArea = float(np.sum(im.array)**2./np.sum(im.array**2.))
427 if image_mask is not None:
428 psfApRadius = max(self.config.minPsfApRadiusPix, 3.0*summary.psfSigma)
429 self.log.debug("Using radius of %.3f (pixels) for psfApFluxDelta metric", psfApRadius)
430 psfTraceRadiusDelta, psfApFluxDelta = compute_psf_image_deltas(
431 image_mask,
432 psf,
433 sampling=self.config.psfGridSampling,
434 ap_radius_pix=psfApRadius,
435 bad_mask_bits=self.config.psfBadMaskPlanes
436 )
437 summary.psfTraceRadiusDelta = float(psfTraceRadiusDelta)
438 summary.psfApFluxDelta = float(psfApFluxDelta)
439 if image_ap_corr_map is not None:
440 if self.config.psfApCorrFieldName not in image_ap_corr_map.keys():
441 self.log.warning(f"{self.config.psfApCorrFieldName} not found in "
442 "image_ap_corr_map. Setting psfApCorrSigmaScaledDelta to NaN.")
443 psfApCorrSigmaScaledDelta = nan
444 else:
445 image_ap_corr_field = image_ap_corr_map[self.config.psfApCorrFieldName]
446 psfApCorrSigmaScaledDelta = compute_ap_corr_sigma_scaled_delta(
447 image_mask,
448 image_ap_corr_field,
449 summary.psfSigma,
450 sampling=self.config.psfGridSampling,
451 bad_mask_bits=self.config.psfBadMaskPlanes,
452 )
453 summary.psfApCorrSigmaScaledDelta = float(psfApCorrSigmaScaledDelta)
455 if sources is None:
456 # No sources are available (as in some tests and rare cases where
457 # the selection criteria in finalizeCharacterization lead to no
458 # good sources).
459 return
461 # Count the total number of psf stars used (prior to stats selection).
462 nPsfStar = sources[self.config.starSelection].sum()
463 summary.nPsfStar = int(nPsfStar)
465 psf_mask = self.starSelector.run(sources).selected
466 nPsfStarsUsedInStats = psf_mask.sum()
468 if nPsfStarsUsedInStats == 0:
469 # No stars to measure statistics, so we must return the defaults
470 # of 0 stars and NaN values.
471 return
473 if sources_is_astropy:
474 psf_cat = sources[psf_mask]
475 else:
476 psf_cat = sources[psf_mask].copy(deep=True)
478 starXX = psf_cat[self.config.starShape + "_xx"]
479 starYY = psf_cat[self.config.starShape + "_yy"]
480 starXY = psf_cat[self.config.starShape + "_xy"]
481 psfXX = psf_cat[self.config.psfShape + "_xx"]
482 psfYY = psf_cat[self.config.psfShape + "_yy"]
483 psfXY = psf_cat[self.config.psfShape + "_xy"]
485 # Use the trace radius for the star size.
486 starSize = np.sqrt(starXX/2. + starYY/2.)
488 starE1 = (starXX - starYY)/(starXX + starYY)
489 starE2 = 2*starXY/(starXX + starYY)
490 starSizeMedian = np.median(starSize)
492 starE = np.sqrt(starE1**2.0 + starE2**2.0)
493 starUnNormalizedE = np.sqrt((starXX - starYY)**2.0 + (2.0*starXY)**2.0)
494 starEMedian = np.median(starE)
495 starUnNormalizedEMedian = np.median(starUnNormalizedE)
496 summary.starEMedian = float(starEMedian)
497 summary.starUnNormalizedEMedian = float(starUnNormalizedEMedian)
499 # Use the trace radius for the psf size.
500 psfSize = np.sqrt(psfXX/2. + psfYY/2.)
501 psfE1 = (psfXX - psfYY)/(psfXX + psfYY)
502 psfE2 = 2*psfXY/(psfXX + psfYY)
504 psfStarDeltaE1Median = np.median(starE1 - psfE1)
505 psfStarDeltaE1Scatter = sigmaMad(starE1 - psfE1, scale="normal")
506 psfStarDeltaE2Median = np.median(starE2 - psfE2)
507 psfStarDeltaE2Scatter = sigmaMad(starE2 - psfE2, scale="normal")
509 psfStarDeltaSizeMedian = np.median(starSize - psfSize)
510 psfStarDeltaSizeScatter = sigmaMad(starSize - psfSize, scale="normal")
511 psfStarScaledDeltaSizeScatter = psfStarDeltaSizeScatter/starSizeMedian
513 summary.psfStarDeltaE1Median = float(psfStarDeltaE1Median)
514 summary.psfStarDeltaE2Median = float(psfStarDeltaE2Median)
515 summary.psfStarDeltaE1Scatter = float(psfStarDeltaE1Scatter)
516 summary.psfStarDeltaE2Scatter = float(psfStarDeltaE2Scatter)
517 summary.psfStarDeltaSizeMedian = float(psfStarDeltaSizeMedian)
518 summary.psfStarDeltaSizeScatter = float(psfStarDeltaSizeScatter)
519 summary.psfStarScaledDeltaSizeScatter = float(psfStarScaledDeltaSizeScatter)
521 # Higher-order moments and derived metrics.
522 if sources_is_astropy:
523 columnNames = psf_cat.colnames
524 else:
525 columnNames = psf_cat.schema.getNames()
526 if self.config.starHigherOrderMomentBase + "_03" in columnNames:
527 starM03 = psf_cat[self.config.starHigherOrderMomentBase + "_03"]
528 starM12 = psf_cat[self.config.starHigherOrderMomentBase + "_12"]
529 starM21 = psf_cat[self.config.starHigherOrderMomentBase + "_21"]
530 starM30 = psf_cat[self.config.starHigherOrderMomentBase + "_30"]
531 starM04 = psf_cat[self.config.starHigherOrderMomentBase + "_04"]
532 starM13 = psf_cat[self.config.starHigherOrderMomentBase + "_13"]
533 starM22 = psf_cat[self.config.starHigherOrderMomentBase + "_22"]
534 starM31 = psf_cat[self.config.starHigherOrderMomentBase + "_31"]
535 starM40 = psf_cat[self.config.starHigherOrderMomentBase + "_40"]
537 starComa1Median = np.nanmedian(starM30 + starM12)
538 starComa2Median = np.nanmedian(starM21 + starM03)
539 starTrefoil1Median = np.nanmedian(starM30 - 3.0*starM12)
540 starTrefoil2Median = np.nanmedian(3.0*starM21 - starM03)
541 starKurtosisMedian = np.nanmedian(starM40 + 2.0*starM22 + starM04)
542 starE41Median = np.nanmedian(starM40 - starM04)
543 starE42Median = np.nanmedian(2.0*(starM31 + starM13))
545 summary.starComa1Median = float(starComa1Median)
546 summary.starComa2Median = float(starComa2Median)
547 summary.starTrefoil1Median = float(starTrefoil1Median)
548 summary.starTrefoil2Median = float(starTrefoil2Median)
549 summary.starKurtosisMedian = float(starKurtosisMedian)
550 summary.starE41Median = float(starE41Median)
551 summary.starE42Median = float(starE42Median)
552 else:
553 self.log.warning("Higher-order moments with base column name %s not found in source "
554 "catalog. Setting the derived metrics (i.e. coma1[2], trefoil1[2], "
555 "kurtosis, e4_1, and e4_2) to nan.", self.config.starHigherOrderMomentBase)
557 if image_mask is not None:
558 maxDistToNearestPsf = maximum_nearest_psf_distance(
559 image_mask,
560 psf_cat,
561 sampling=self.config.psfSampling,
562 bad_mask_bits=self.config.psfBadMaskPlanes
563 )
564 summary.maxDistToNearestPsf = float(maxDistToNearestPsf)
566 def update_wcs_stats(self, summary, wcs, bbox, visitInfo):
567 """Compute all summary-statistic fields that depend on the WCS model.
569 Parameters
570 ----------
571 summary : `lsst.afw.image.ExposureSummaryStats`
572 Summary object to update in-place.
573 wcs : `lsst.afw.geom.SkyWcs` or `None`
574 Astrometric calibration model. If `None`, all fields that depend
575 on the WCS will be reset (generally to NaN).
576 bbox : `lsst.geom.Box2I`
577 Bounding box of the image for which summary stats are being
578 computed.
579 visitInfo : `lsst.afw.image.VisitInfo`
580 Observation information used in together with ``wcs`` to compute
581 the zenith distance.
582 """
583 nan = float("nan")
584 summary.raCorners = [nan]*4
585 summary.decCorners = [nan]*4
586 summary.ra = nan
587 summary.dec = nan
588 summary.pixelScale = nan
589 summary.zenithDistance = nan
591 if wcs is None:
592 return
594 sph_pts = wcs.pixelToSky(geom.Box2D(bbox).getCorners())
595 summary.raCorners = [float(sph.getRa().asDegrees()) for sph in sph_pts]
596 summary.decCorners = [float(sph.getDec().asDegrees()) for sph in sph_pts]
598 sph_pt = wcs.pixelToSky(bbox.getCenter())
599 summary.ra = sph_pt.getRa().asDegrees()
600 summary.dec = sph_pt.getDec().asDegrees()
601 summary.pixelScale = wcs.getPixelScale(bbox.getCenter()).asArcseconds()
603 date = visitInfo.getDate()
605 if date.isValid():
606 # We compute the zenithDistance at the center of the detector
607 # rather than use the boresight value available via the visitInfo,
608 # because the zenithDistance may vary significantly over a large
609 # field of view.
610 observatory = visitInfo.getObservatory()
611 loc = EarthLocation(lat=observatory.getLatitude().asDegrees()*units.deg,
612 lon=observatory.getLongitude().asDegrees()*units.deg,
613 height=observatory.getElevation()*units.m)
614 obstime = Time(visitInfo.getDate().get(system=DateTime.MJD),
615 location=loc, format="mjd")
616 coord = SkyCoord(
617 summary.ra*units.degree,
618 summary.dec*units.degree,
619 obstime=obstime,
620 location=loc,
621 )
622 with warnings.catch_warnings():
623 warnings.simplefilter("ignore")
624 altaz = coord.transform_to(AltAz)
626 summary.zenithDistance = float(90.0 - altaz.alt.degree)
628 def update_photo_calib_stats(self, summary, photo_calib):
629 """Compute all summary-statistic fields that depend on the photometric
630 calibration model.
632 Parameters
633 ----------
634 summary : `lsst.afw.image.ExposureSummaryStats`
635 Summary object to update in-place.
636 photo_calib : `lsst.afw.image.PhotoCalib` or `None`
637 Photometric calibration model. If `None`, all fields that depend
638 on the photometric calibration will be reset (generally to NaN).
639 """
640 if photo_calib is not None:
641 summary.zeroPoint = float(2.5*np.log10(photo_calib.getInstFluxAtZeroMagnitude()))
642 else:
643 summary.zeroPoint = float("nan")
645 def update_background_stats(self, summary, background):
646 """Compute summary-statistic fields that depend only on the
647 background model.
649 Parameters
650 ----------
651 summary : `lsst.afw.image.ExposureSummaryStats`
652 Summary object to update in-place.
653 background : `lsst.afw.math.BackgroundList` or `None`
654 Background model. If `None`, all fields that depend on the
655 background will be reset (generally to NaN).
657 Notes
658 -----
659 This does not include fields that depend on the background-subtracted
660 masked image; when the background changes, it should generally be
661 applied to the image and `update_masked_image_stats` should be called
662 as well.
663 """
664 if background is not None:
665 bgStats = (bg[0].getStatsImage().getImage().array
666 for bg in background)
667 summary.skyBg = float(sum(np.median(bg[np.isfinite(bg)]) for bg in bgStats))
668 else:
669 summary.skyBg = float("nan")
671 def update_masked_image_stats(self, summary, masked_image):
672 """Compute summary-statistic fields that depend on the masked image
673 itself.
675 Parameters
676 ----------
677 summary : `lsst.afw.image.ExposureSummaryStats`
678 Summary object to update in-place.
679 masked_image : `lsst.afw.image.MaskedImage` or `None`
680 Masked image. If `None`, all fields that depend
681 on the masked image will be reset (generally to NaN).
682 """
683 nan = float("nan")
684 if masked_image is None:
685 summary.skyNoise = nan
686 summary.meanVar = nan
687 return
688 statsCtrl = afwMath.StatisticsControl()
689 statsCtrl.setNumSigmaClip(self.config.sigmaClip)
690 statsCtrl.setNumIter(self.config.clipIter)
691 statsCtrl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.badMaskPlanes))
692 statsCtrl.setNanSafe(True)
694 statObj = afwMath.makeStatistics(masked_image, afwMath.STDEVCLIP, statsCtrl)
695 skyNoise, _ = statObj.getResult(afwMath.STDEVCLIP)
696 summary.skyNoise = skyNoise
698 statObj = afwMath.makeStatistics(masked_image.variance, masked_image.mask, afwMath.MEANCLIP,
699 statsCtrl)
700 meanVar, _ = statObj.getResult(afwMath.MEANCLIP)
701 summary.meanVar = meanVar
703 def update_effective_time_stats(self, summary, exposure):
704 """Compute effective exposure time statistics to estimate depth.
706 The effective exposure time is the equivalent shutter open
707 time that would be needed under nominal conditions to give the
708 same signal-to-noise for a point source as what is achieved by
709 the observation of interest. This metric combines measurements
710 of the point-spread function, the sky brightness, and the
711 transparency. It assumes that the observation is
712 sky-background dominated.
714 .. _teff_definitions:
716 The effective exposure time and its subcomponents are defined in [1]_.
718 References
719 ----------
721 .. [1] Neilsen, E.H., Bernstein, G., Gruendl, R., and Kent, S. (2016).
722 Limiting Magnitude, \tau, teff, and Image Quality in DES Year 1
723 https://www.osti.gov/biblio/1250877/
725 Parameters
726 ----------
727 summary : `lsst.afw.image.ExposureSummaryStats`
728 Summary object to update in-place.
729 exposure : `lsst.afw.image.ExposureF`
730 Exposure to grab band and exposure time metadata
732 """
733 nan = float("nan")
734 summary.effTime = nan
735 summary.effTimePsfSigmaScale = nan
736 summary.effTimeSkyBgScale = nan
737 summary.effTimeZeroPointScale = nan
739 exposureTime = exposure.getInfo().getVisitInfo().getExposureTime()
740 filterLabel = exposure.getFilter()
741 if (filterLabel is None) or (not filterLabel.hasBandLabel):
742 band = None
743 else:
744 band = filterLabel.bandLabel
746 if band is None:
747 self.log.warning("No band associated with exposure; effTime not calculated.")
748 return
750 # PSF component
751 if np.isnan(summary.psfSigma):
752 self.log.debug("PSF sigma is NaN")
753 f_eff = nan
754 elif band not in self.config.fiducialPsfSigma:
755 self.log.debug(f"Fiducial PSF value not found for {band}")
756 f_eff = nan
757 else:
758 fiducialPsfSigma = self.config.fiducialPsfSigma[band]
759 f_eff = (summary.psfSigma / fiducialPsfSigma)**-2
761 # Transparency component
762 # Note: Assumes that the zeropoint includes the exposure time
763 if np.isnan(summary.zeroPoint):
764 self.log.debug("Zero point is NaN")
765 c_eff = nan
766 elif band not in self.config.fiducialZeroPoint:
767 self.log.debug(f"Fiducial zero point value not found for {band}")
768 c_eff = nan
769 else:
770 fiducialZeroPoint = self.config.fiducialZeroPoint[band]
771 zeroPointDiff = fiducialZeroPoint - (summary.zeroPoint - 2.5*np.log10(exposureTime))
772 c_eff = min(10**(-2.0*(zeroPointDiff)/2.5), self.config.maxEffectiveTransparency)
774 # Sky brightness component (convert to cts/s)
775 if np.isnan(summary.skyBg):
776 self.log.debug("Sky background is NaN")
777 b_eff = nan
778 elif band not in self.config.fiducialSkyBackground:
779 self.log.debug(f"Fiducial sky background value not found for {band}")
780 b_eff = nan
781 else:
782 fiducialSkyBackground = self.config.fiducialSkyBackground[band]
783 b_eff = fiducialSkyBackground/(summary.skyBg/exposureTime)
785 # Old effective exposure time scale factor
786 # t_eff = f_eff * c_eff * b_eff
788 # New effective exposure time (seconds) from magnitude limit
789 if band not in self.config.fiducialMagLim:
790 self.log.debug(f"Fiducial magnitude limit not found for {band}")
791 effectiveTime = nan
792 elif band not in self.config.fiducialExpTime:
793 self.log.debug(f"Fiducial exposure time not found for {band}")
794 effectiveTime = nan
795 else:
796 effectiveTime = compute_effective_time(summary.magLim,
797 self.config.fiducialMagLim[band],
798 self.config.fiducialExpTime[band])
800 # Output quantities
801 summary.effTime = float(effectiveTime)
802 summary.effTimePsfSigmaScale = float(f_eff)
803 summary.effTimeSkyBgScale = float(b_eff)
804 summary.effTimeZeroPointScale = float(c_eff)
806 def update_magnitude_limit_stats(self, summary, exposure):
807 """Compute a summary statistic for the point-source magnitude limit at
808 a given signal-to-noise ratio using exposure-level metadata.
810 The magnitude limit is calculated at a given SNR from the PSF,
811 sky, zeropoint, and readnoise in accordance with SMTN-002 [1]_,
812 LSE-40 [2]_, and DMTN-296 [3]_.
814 References
815 ----------
817 .. [1] "Calculating LSST limiting magnitudes and SNR" (SMTN-002)
818 .. [2] "Photon Rates and SNR Calculations" (LSE-40)
819 .. [3] "Calculations of Image and Catalog Depth" (DMTN-296)
822 Parameters
823 ----------
824 summary : `lsst.afw.image.ExposureSummaryStats`
825 Summary object to update in-place.
826 exposure : `lsst.afw.image.ExposureF`
827 Exposure to grab band and exposure time metadata
828 """
829 if exposure.getDetector() is None:
830 summary.magLim = float("nan")
831 return
833 # Calculate the average readnoise [e-]
834 readNoiseList = list(ipIsr.getExposureReadNoises(exposure).values())
835 readNoise = np.nanmean(readNoiseList)
836 if np.isnan(readNoise):
837 readNoise = 0.0
838 self.log.warning("Read noise set to NaN! Setting readNoise to 0.0 to proceed.")
840 # Calculate the average gain [e-/ADU]
841 gainList = list(ipIsr.getExposureGains(exposure).values())
842 gain = np.nanmean(gainList)
843 if np.isnan(gain):
844 self.log.warning("Gain set to NaN! Setting magLim to NaN.")
845 summary.magLim = float("nan")
846 return
848 # Get the image units (default to "adu" if metadata key absent)
849 md = exposure.getMetadata()
850 if md.get("LSST ISR UNITS", "adu") == "electron":
851 gain = 1.0
853 # Convert readNoise to image units
854 readNoise /= gain
856 # Calculate the limiting magnitude.
857 # Note 1: Assumes that the image and readnoise have the same units
858 # Note 2: Assumes that the zeropoint includes the exposure time
859 magLim = compute_magnitude_limit(summary.psfArea, summary.skyBg,
860 summary.zeroPoint, readNoise,
861 gain, self.config.magLimSnr)
863 summary.magLim = float(magLim)
866def maximum_nearest_psf_distance(
867 image_mask,
868 psf_cat,
869 sampling=8,
870 bad_mask_bits=["BAD", "CR", "INTRP", "SAT", "SUSPECT", "NO_DATA", "EDGE"],
871):
872 """Compute the maximum distance of an unmasked pixel to its nearest PSF.
874 Parameters
875 ----------
876 image_mask : `lsst.afw.image.Mask`
877 The mask plane associated with the exposure.
878 psf_cat : `lsst.afw.table.SourceCatalog` or `astropy.table.Table`
879 Catalog containing only the stars used in the PSF modeling.
880 sampling : `int`
881 Sampling rate in each dimension to create the grid of points on which
882 to evaluate the distance to the nearest PSF star. The tradeoff is
883 between adequate sampling versus speed.
884 bad_mask_bits : `list` [`str`]
885 Mask bits required to be absent for a pixel to be considered
886 "unmasked".
888 Returns
889 -------
890 max_dist_to_nearest_psf : `float`
891 The maximum distance (in pixels) of an unmasked pixel to its nearest
892 PSF model star.
893 """
894 mask_arr = image_mask.array[::sampling, ::sampling]
895 bitmask = image_mask.getPlaneBitMask(bad_mask_bits)
896 good = ((mask_arr & bitmask) == 0)
898 x = np.arange(good.shape[1]) * sampling
899 y = np.arange(good.shape[0]) * sampling
900 xx, yy = np.meshgrid(x, y)
902 dist_to_nearest_psf = np.full(good.shape, np.inf)
903 for psf in psf_cat:
904 x_psf = psf["slot_Centroid_x"]
905 y_psf = psf["slot_Centroid_y"]
906 dist_to_nearest_psf = np.minimum(dist_to_nearest_psf, np.hypot(xx - x_psf, yy - y_psf))
907 unmasked_dists = dist_to_nearest_psf * good
908 max_dist_to_nearest_psf = np.max(unmasked_dists)
910 return max_dist_to_nearest_psf
913def compute_psf_image_deltas(
914 image_mask,
915 image_psf,
916 sampling=96,
917 ap_radius_pix=3.0,
918 bad_mask_bits=["BAD", "CR", "INTRP", "SAT", "SUSPECT", "NO_DATA", "EDGE"],
919):
920 """Compute the delta between the maximum and minimum model PSF trace radius
921 values evaluated on a grid of points lying in the unmasked region of the
922 image.
924 Parameters
925 ----------
926 image_mask : `lsst.afw.image.Mask`
927 The mask plane associated with the exposure.
928 image_psf : `lsst.afw.detection.Psf`
929 The PSF model associated with the exposure.
930 sampling : `int`, optional
931 Sampling rate in each dimension to create the grid of points at which
932 to evaluate ``image_psf``s trace radius value. The tradeoff is between
933 adequate sampling versus speed.
934 ap_radius_pix : `float`, optional
935 Radius in pixels of the aperture on which to measure the flux of the
936 PSF model.
937 bad_mask_bits : `list` [`str`], optional
938 Mask bits required to be absent for a pixel to be considered
939 "unmasked".
941 Returns
942 -------
943 psf_trace_radius_delta, psf_ap_flux_delta : `float`
944 The delta (in pixels) between the maximum and minimum model PSF trace
945 radius values and the PSF aperture fluxes (with aperture radius of
946 max(2, 3*psfSigma)) evaluated on the x,y-grid subsampled on the
947 unmasked detector pixels by a factor of ``sampling``. If both the
948 model PSF trace radius value and aperture flux value on the grid
949 evaluate to NaN, then NaNs are returned immediately.
950 """
951 psf_trace_radius_list = []
952 psf_ap_flux_list = []
953 mask_arr = image_mask.array[::sampling, ::sampling]
954 bitmask = image_mask.getPlaneBitMask(bad_mask_bits)
955 good = ((mask_arr & bitmask) == 0)
957 x = np.arange(good.shape[1]) * sampling
958 y = np.arange(good.shape[0]) * sampling
959 xx, yy = np.meshgrid(x, y)
961 for x_mesh, y_mesh, good_mesh in zip(xx, yy, good):
962 for x_point, y_point, is_good in zip(x_mesh, y_mesh, good_mesh):
963 if is_good:
964 position = geom.Point2D(x_point, y_point)
965 psf_trace_radius = image_psf.computeShape(position).getTraceRadius()
966 psf_ap_flux = image_psf.computeApertureFlux(ap_radius_pix, position)
967 if ~np.isfinite(psf_trace_radius) and ~np.isfinite(psf_ap_flux):
968 return float("nan"), float("nan")
969 psf_trace_radius_list.append(psf_trace_radius)
970 psf_ap_flux_list.append(psf_ap_flux)
972 psf_trace_radius_delta = np.max(psf_trace_radius_list) - np.min(psf_trace_radius_list)
973 if np.any(np.asarray(psf_ap_flux_list) < 0.0): # Consider any -ve flux entry as "bad".
974 psf_ap_flux_delta = float("nan")
975 else:
976 psf_ap_flux_delta = np.max(psf_ap_flux_list) - np.min(psf_ap_flux_list)
978 return psf_trace_radius_delta, psf_ap_flux_delta
981def compute_ap_corr_sigma_scaled_delta(
982 image_mask,
983 image_ap_corr_field,
984 psfSigma,
985 sampling=96,
986 bad_mask_bits=["BAD", "CR", "INTRP", "SAT", "SUSPECT", "NO_DATA", "EDGE"],
987):
988 """Compute the delta between the maximum and minimum aperture correction
989 values scaled (divided) by ``psfSigma`` for the given field representation,
990 ``image_ap_corr_field`` evaluated on a grid of points lying in the
991 unmasked region of the image.
993 Parameters
994 ----------
995 image_mask : `lsst.afw.image.Mask`
996 The mask plane associated with the exposure.
997 image_ap_corr_field : `lsst.afw.math.ChebyshevBoundedField`
998 The ChebyshevBoundedField representation of the aperture correction
999 of interest for the exposure.
1000 psfSigma : `float`
1001 The PSF model second-moments determinant radius (center of chip)
1002 in pixels.
1003 sampling : `int`, optional
1004 Sampling rate in each dimension to create the grid of points at which
1005 to evaluate ``image_psf``s trace radius value. The tradeoff is between
1006 adequate sampling versus speed.
1007 bad_mask_bits : `list` [`str`], optional
1008 Mask bits required to be absent for a pixel to be considered
1009 "unmasked".
1011 Returns
1012 -------
1013 ap_corr_sigma_scaled_delta : `float`
1014 The delta between the maximum and minimum of the (multiplicative)
1015 aperture correction values scaled (divided) by ``psfSigma`` evaluated
1016 on the x,y-grid subsampled on the unmasked detector pixels by a factor
1017 of ``sampling``. If the aperture correction evaluates to NaN on any
1018 of the grid points, this is set to NaN.
1019 """
1020 mask_arr = image_mask.array[::sampling, ::sampling]
1021 bitmask = image_mask.getPlaneBitMask(bad_mask_bits)
1022 good = ((mask_arr & bitmask) == 0)
1024 x = np.arange(good.shape[1], dtype=np.float64) * sampling
1025 y = np.arange(good.shape[0], dtype=np.float64) * sampling
1026 xx, yy = np.meshgrid(x, y)
1028 ap_corr = image_ap_corr_field.evaluate(xx.ravel(), yy.ravel()).reshape(xx.shape)
1029 ap_corr_good = ap_corr[good]
1030 if ~np.isfinite(ap_corr_good).all():
1031 ap_corr_sigma_scaled_delta = float("nan")
1032 else:
1033 ap_corr_sigma_scaled_delta = (np.max(ap_corr_good) - np.min(ap_corr_good))/psfSigma
1035 return ap_corr_sigma_scaled_delta
1038def compute_magnitude_limit(
1039 psfArea,
1040 skyBg,
1041 zeroPoint,
1042 readNoise,
1043 gain,
1044 snr
1045):
1046 """Compute the expected point-source magnitude limit at a given
1047 signal-to-noise ratio given the exposure-level metadata. Based on
1048 the signal-to-noise formula provided in SMTN-002 (see LSE-40 for
1049 more details on the calculation).
1051 SNR = C / sqrt( C/g + (B/g + sigma_inst**2) * neff )
1053 where C is the counts from the source, B is counts from the (sky)
1054 background, sigma_inst is the instrumental (read) noise, neff is
1055 the effective size of the PSF, and g is the gain in e-/ADU. Note
1056 that input values of ``skyBg``, ``zeroPoint``, and ``readNoise``
1057 should all consistently be in electrons or ADU.
1059 Parameters
1060 ----------
1061 psfArea : `float`
1062 The effective area of the PSF [pix].
1063 skyBg : `float`
1064 The sky background counts for the exposure [ADU or e-].
1065 zeroPoint : `float`
1066 The zeropoint (includes exposure time) [ADU or e-].
1067 readNoise : `float`
1068 The instrumental read noise for the exposure [ADU or e-].
1069 gain : `float`
1070 The instrumental gain for the exposure [e-/ADU]. The gain should
1071 be 1.0 if the skyBg, zeroPoint, and readNoise are in e-.
1072 snr : `float`
1073 Signal-to-noise ratio at which magnitude limit is calculated.
1075 Returns
1076 -------
1077 magnitude_limit : `float`
1078 The expected magnitude limit at the given signal to noise.
1079 """
1080 # Effective number of pixels within the PSF calculated directly
1081 # from the PSF model.
1082 neff = psfArea
1084 # Instrumental noise (read noise only)
1085 sigma_inst = readNoise
1087 # Total background counts derived from Eq. 45 of LSE-40
1088 background = (skyBg/gain + sigma_inst**2) * neff
1089 # Source counts to achieve the desired SNR (from quadratic formula)
1090 # Note typo in gain exponent in Eq. 45 of LSE-40 (DM-53234)
1091 source = (snr**2)/(2*gain) + np.sqrt((snr**4)/(4*gain**2) + snr**2 * background)
1093 # Convert to a magnitude using the zeropoint.
1094 # Note: Zeropoint currently includes exposure time
1095 magnitude_limit = -2.5*np.log10(source) + zeroPoint
1097 return magnitude_limit
1100def psf_sigma_to_psf_area(psfSigma, pixelScale):
1101 """Convert psfSigma [pixels] to an approximate psfArea [pixel^2] as defined in LSE-40.
1102 This is the same implementation followed by SMTN-002 to estimate SNR=5 magnitude limits.
1104 Parameters
1105 ----------
1106 psfSigma : `float`
1107 The PSF sigma value [pix].
1108 pixelScale : `float`
1109 The pixel scale [arcsec/pix].
1111 Returns
1112 -------
1113 psf_area : `float`
1114 Approximation of the PSF area [pix^2]
1115 """
1116 # Follow SMTN-002 to convert to geometric and effective FWHM
1117 fwhm_geom = psfSigma * 2.355 * pixelScale
1118 fwhm_eff = (fwhm_geom - 0.052) / 0.822
1119 # Area prefactor comes from LSE-40
1120 psf_area = 2.266 * (fwhm_eff / pixelScale)**2
1121 return psf_area
1124def compute_effective_time(
1125 magLim,
1126 fiducialMagLim,
1127 fiducialExpTime
1128):
1129 """Compute the effective exposure time from m5 following the prescription described in SMTN-296.
1131 teff = 10**(0.8 * (magLim - fiducialMagLim) ) * fiducialExpTime
1133 where `magLim` is the magnitude limit, `fiducialMagLim` is the fiducial magnitude limit calculated from
1134 the fiducial values of the ``psfArea``, ``skyBg``, ``zeroPoint``, and ``readNoise``, and
1135 `fiducialExpTime` is the fiducial exposure time (s).
1137 Parameters
1138 ----------
1139 magLim : `float`
1140 The measured magnitude limit [mag].
1141 fiducialMagLim : `float`
1142 The fiducial magnitude limit [mag].
1143 fiducialExpTime : `float`
1144 The fiducial exposure time [s].
1146 Returns
1147 -------
1148 effectiveTime : `float`
1149 The effective exposure time.
1150 """
1151 effectiveTime = 10**(0.8 * (magLim - fiducialMagLim)) * fiducialExpTime
1152 return effectiveTime