Coverage for python / lsst / pipe / tasks / computeExposureSummaryStats.py: 12%
417 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-13 01:56 -0700
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, summary=None):
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`
310 summary : `lsst.afw.image.ExposureSummary`, optional
311 Summary object to be appended to if not `None`. If `None` a new
312 summary object will be created.
314 Returns
315 -------
316 summary : `lsst.afw.image.ExposureSummary`
317 """
318 self.log.info("Measuring exposure statistics")
320 if summary is None:
321 summary = afwImage.ExposureSummaryStats()
323 # Set exposure time.
324 exposureTime = exposure.getInfo().getVisitInfo().getExposureTime()
325 summary.expTime = exposureTime
327 bbox = exposure.getBBox()
329 psf = exposure.getPsf()
330 self.update_psf_stats(
331 summary, psf, bbox, sources, image_mask=exposure.mask, image_ap_corr_map=exposure.apCorrMap
332 )
334 wcs = exposure.getWcs()
335 visitInfo = exposure.getInfo().getVisitInfo()
336 self.update_wcs_stats(summary, wcs, bbox, visitInfo)
338 photoCalib = exposure.getPhotoCalib()
339 self.update_photo_calib_stats(summary, photoCalib)
341 self.update_background_stats(summary, background)
343 self.update_masked_image_stats(summary, exposure.getMaskedImage())
345 self.update_magnitude_limit_stats(summary, exposure)
347 self.update_effective_time_stats(summary, exposure)
349 md = exposure.getMetadata()
350 if "PSF_ADAPTIVE_THRESHOLD_VALUE" in md:
351 summary.psfAdaptiveThresholdValue = md["PSF_ADAPTIVE_THRESHOLD_VALUE"]
352 if "PSF_ADAPTIVE_INCLUDE_THRESHOLD_MULTIPLIER" in md:
353 summary.psfAdaptiveIncludeThresholdMultiplier = md["PSF_ADAPTIVE_INCLUDE_THRESHOLD_MULTIPLIER"]
354 if "SFM_ASTROM_OFFSET_MEAN" in md:
355 summary.astromOffsetMean = md["SFM_ASTROM_OFFSET_MEAN"]
356 summary.astromOffsetStd = md["SFM_ASTROM_OFFSET_STD"]
358 return summary
360 def update_psf_stats(
361 self,
362 summary,
363 psf,
364 bbox,
365 sources=None,
366 image_mask=None,
367 image_ap_corr_map=None,
368 sources_is_astropy=False,
369 ):
370 """Compute all summary-statistic fields that depend on the PSF model.
372 Parameters
373 ----------
374 summary : `lsst.afw.image.ExposureSummaryStats`
375 Summary object to update in-place.
376 psf : `lsst.afw.detection.Psf` or `None`
377 Point spread function model. If `None`, all fields that depend on
378 the PSF will be reset (generally to NaN).
379 bbox : `lsst.geom.Box2I`
380 Bounding box of the image for which summary stats are being
381 computed.
382 sources : `lsst.afw.table.SourceCatalog` or `astropy.table.Table`
383 Catalog for quantities that are computed from source table columns.
384 If `None`, these quantities will be reset (generally to NaN).
385 The type of this table must correspond to the
386 ``sources_is_astropy`` argument.
387 image_mask : `lsst.afw.image.Mask`, optional
388 Mask image that may be used to compute distance-to-nearest-star
389 metrics.
390 sources_is_astropy : `bool`, optional
391 Whether ``sources`` is an `astropy.table.Table` instance instead
392 of an `lsst.afw.table.Catalog` instance. Default is `False` (the
393 latter).
394 """
395 nan = float("nan")
396 summary.psfSigma = nan
397 summary.psfIxx = nan
398 summary.psfIyy = nan
399 summary.psfIxy = nan
400 summary.psfArea = nan
401 summary.nPsfStar = 0
402 summary.psfStarDeltaE1Median = nan
403 summary.psfStarDeltaE2Median = nan
404 summary.psfStarDeltaE1Scatter = nan
405 summary.psfStarDeltaE2Scatter = nan
406 summary.psfStarDeltaSizeMedian = nan
407 summary.psfStarDeltaSizeScatter = nan
408 summary.psfStarScaledDeltaSizeScatter = nan
409 summary.maxDistToNearestPsf = nan
410 summary.psfTraceRadiusDelta = nan
411 summary.psfApFluxDelta = nan
412 summary.psfApCorrSigmaScaledDelta = nan
413 summary.starEMedian = nan
414 summary.starUnNormalizedEMedian = nan
415 summary.starComa1Median = nan
416 summary.starComa2Median = nan
417 summary.starTrefoil1Median = nan
418 summary.starTrefoil2Median = nan
419 summary.starKurtosisMedian = nan
420 summary.starE41Median = nan
421 summary.starE42Median = nan
423 if psf is None:
424 return
425 shape = psf.computeShape(bbox.getCenter())
426 summary.psfSigma = shape.getDeterminantRadius()
427 summary.psfIxx = shape.getIxx()
428 summary.psfIyy = shape.getIyy()
429 summary.psfIxy = shape.getIxy()
430 im = psf.computeKernelImage(bbox.getCenter())
431 # The calculation of effective psf area is taken from
432 # ls.st/srd Equation 1.
433 summary.psfArea = float(np.sum(im.array)**2./np.sum(im.array**2.))
435 if image_mask is not None:
436 psfApRadius = max(self.config.minPsfApRadiusPix, 3.0*summary.psfSigma)
437 self.log.debug("Using radius of %.3f (pixels) for psfApFluxDelta metric", psfApRadius)
438 psfTraceRadiusDelta, psfApFluxDelta = compute_psf_image_deltas(
439 image_mask,
440 psf,
441 sampling=self.config.psfGridSampling,
442 ap_radius_pix=psfApRadius,
443 bad_mask_bits=self.config.psfBadMaskPlanes
444 )
445 summary.psfTraceRadiusDelta = float(psfTraceRadiusDelta)
446 summary.psfApFluxDelta = float(psfApFluxDelta)
447 if image_ap_corr_map is not None:
448 if self.config.psfApCorrFieldName not in image_ap_corr_map.keys():
449 self.log.warning(f"{self.config.psfApCorrFieldName} not found in "
450 "image_ap_corr_map. Setting psfApCorrSigmaScaledDelta to NaN.")
451 psfApCorrSigmaScaledDelta = nan
452 else:
453 image_ap_corr_field = image_ap_corr_map[self.config.psfApCorrFieldName]
454 psfApCorrSigmaScaledDelta = compute_ap_corr_sigma_scaled_delta(
455 image_mask,
456 image_ap_corr_field,
457 summary.psfSigma,
458 sampling=self.config.psfGridSampling,
459 bad_mask_bits=self.config.psfBadMaskPlanes,
460 )
461 summary.psfApCorrSigmaScaledDelta = float(psfApCorrSigmaScaledDelta)
463 if sources is None:
464 # No sources are available (as in some tests and rare cases where
465 # the selection criteria in finalizeCharacterization lead to no
466 # good sources).
467 return
469 # Count the total number of psf stars used (prior to stats selection).
470 nPsfStar = sources[self.config.starSelection].sum()
471 summary.nPsfStar = int(nPsfStar)
473 psf_mask = self.starSelector.run(sources).selected
474 nPsfStarsUsedInStats = psf_mask.sum()
476 if nPsfStarsUsedInStats == 0:
477 # No stars to measure statistics, so we must return the defaults
478 # of 0 stars and NaN values.
479 return
481 if sources_is_astropy:
482 psf_cat = sources[psf_mask]
483 else:
484 psf_cat = sources[psf_mask].copy(deep=True)
486 starXX = psf_cat[self.config.starShape + "_xx"]
487 starYY = psf_cat[self.config.starShape + "_yy"]
488 starXY = psf_cat[self.config.starShape + "_xy"]
489 psfXX = psf_cat[self.config.psfShape + "_xx"]
490 psfYY = psf_cat[self.config.psfShape + "_yy"]
491 psfXY = psf_cat[self.config.psfShape + "_xy"]
493 # Use the trace radius for the star size.
494 starSize = np.sqrt(starXX/2. + starYY/2.)
496 starE1 = (starXX - starYY)/(starXX + starYY)
497 starE2 = 2*starXY/(starXX + starYY)
498 starSizeMedian = np.median(starSize)
500 starE = np.sqrt(starE1**2.0 + starE2**2.0)
501 starUnNormalizedE = np.sqrt((starXX - starYY)**2.0 + (2.0*starXY)**2.0)
502 starEMedian = np.median(starE)
503 starUnNormalizedEMedian = np.median(starUnNormalizedE)
504 summary.starEMedian = float(starEMedian)
505 summary.starUnNormalizedEMedian = float(starUnNormalizedEMedian)
507 # Use the trace radius for the psf size.
508 psfSize = np.sqrt(psfXX/2. + psfYY/2.)
509 psfE1 = (psfXX - psfYY)/(psfXX + psfYY)
510 psfE2 = 2*psfXY/(psfXX + psfYY)
512 psfStarDeltaE1Median = np.median(starE1 - psfE1)
513 psfStarDeltaE1Scatter = sigmaMad(starE1 - psfE1, scale="normal")
514 psfStarDeltaE2Median = np.median(starE2 - psfE2)
515 psfStarDeltaE2Scatter = sigmaMad(starE2 - psfE2, scale="normal")
517 psfStarDeltaSizeMedian = np.median(starSize - psfSize)
518 psfStarDeltaSizeScatter = sigmaMad(starSize - psfSize, scale="normal")
519 psfStarScaledDeltaSizeScatter = psfStarDeltaSizeScatter/starSizeMedian
521 summary.psfStarDeltaE1Median = float(psfStarDeltaE1Median)
522 summary.psfStarDeltaE2Median = float(psfStarDeltaE2Median)
523 summary.psfStarDeltaE1Scatter = float(psfStarDeltaE1Scatter)
524 summary.psfStarDeltaE2Scatter = float(psfStarDeltaE2Scatter)
525 summary.psfStarDeltaSizeMedian = float(psfStarDeltaSizeMedian)
526 summary.psfStarDeltaSizeScatter = float(psfStarDeltaSizeScatter)
527 summary.psfStarScaledDeltaSizeScatter = float(psfStarScaledDeltaSizeScatter)
529 # Higher-order moments and derived metrics.
530 if sources_is_astropy:
531 columnNames = psf_cat.colnames
532 else:
533 columnNames = psf_cat.schema.getNames()
534 if self.config.starHigherOrderMomentBase + "_03" in columnNames:
535 starM03 = psf_cat[self.config.starHigherOrderMomentBase + "_03"]
536 starM12 = psf_cat[self.config.starHigherOrderMomentBase + "_12"]
537 starM21 = psf_cat[self.config.starHigherOrderMomentBase + "_21"]
538 starM30 = psf_cat[self.config.starHigherOrderMomentBase + "_30"]
539 starM04 = psf_cat[self.config.starHigherOrderMomentBase + "_04"]
540 starM13 = psf_cat[self.config.starHigherOrderMomentBase + "_13"]
541 starM22 = psf_cat[self.config.starHigherOrderMomentBase + "_22"]
542 starM31 = psf_cat[self.config.starHigherOrderMomentBase + "_31"]
543 starM40 = psf_cat[self.config.starHigherOrderMomentBase + "_40"]
545 starComa1Median = np.nanmedian(starM30 + starM12)
546 starComa2Median = np.nanmedian(starM21 + starM03)
547 starTrefoil1Median = np.nanmedian(starM30 - 3.0*starM12)
548 starTrefoil2Median = np.nanmedian(3.0*starM21 - starM03)
549 starKurtosisMedian = np.nanmedian(starM40 + 2.0*starM22 + starM04)
550 starE41Median = np.nanmedian(starM40 - starM04)
551 starE42Median = np.nanmedian(2.0*(starM31 + starM13))
553 summary.starComa1Median = float(starComa1Median)
554 summary.starComa2Median = float(starComa2Median)
555 summary.starTrefoil1Median = float(starTrefoil1Median)
556 summary.starTrefoil2Median = float(starTrefoil2Median)
557 summary.starKurtosisMedian = float(starKurtosisMedian)
558 summary.starE41Median = float(starE41Median)
559 summary.starE42Median = float(starE42Median)
560 else:
561 self.log.warning("Higher-order moments with base column name %s not found in source "
562 "catalog. Setting the derived metrics (i.e. coma1[2], trefoil1[2], "
563 "kurtosis, e4_1, and e4_2) to nan.", self.config.starHigherOrderMomentBase)
565 if image_mask is not None:
566 maxDistToNearestPsf = maximum_nearest_psf_distance(
567 image_mask,
568 psf_cat,
569 sampling=self.config.psfSampling,
570 bad_mask_bits=self.config.psfBadMaskPlanes
571 )
572 summary.maxDistToNearestPsf = float(maxDistToNearestPsf)
574 def update_wcs_stats(self, summary, wcs, bbox, visitInfo):
575 """Compute all summary-statistic fields that depend on the WCS model.
577 Parameters
578 ----------
579 summary : `lsst.afw.image.ExposureSummaryStats`
580 Summary object to update in-place.
581 wcs : `lsst.afw.geom.SkyWcs` or `None`
582 Astrometric calibration model. If `None`, all fields that depend
583 on the WCS will be reset (generally to NaN).
584 bbox : `lsst.geom.Box2I`
585 Bounding box of the image for which summary stats are being
586 computed.
587 visitInfo : `lsst.afw.image.VisitInfo`
588 Observation information used in together with ``wcs`` to compute
589 the zenith distance.
590 """
591 nan = float("nan")
592 summary.raCorners = [nan]*4
593 summary.decCorners = [nan]*4
594 summary.ra = nan
595 summary.dec = nan
596 summary.pixelScale = nan
597 summary.zenithDistance = nan
599 if wcs is None:
600 return
602 sph_pts = wcs.pixelToSky(geom.Box2D(bbox).getCorners())
603 summary.raCorners = [float(sph.getRa().asDegrees()) for sph in sph_pts]
604 summary.decCorners = [float(sph.getDec().asDegrees()) for sph in sph_pts]
606 sph_pt = wcs.pixelToSky(bbox.getCenter())
607 summary.ra = sph_pt.getRa().asDegrees()
608 summary.dec = sph_pt.getDec().asDegrees()
609 summary.pixelScale = wcs.getPixelScale(bbox.getCenter()).asArcseconds()
611 date = visitInfo.getDate()
613 if date.isValid():
614 # We compute the zenithDistance at the center of the detector
615 # rather than use the boresight value available via the visitInfo,
616 # because the zenithDistance may vary significantly over a large
617 # field of view.
618 observatory = visitInfo.getObservatory()
619 loc = EarthLocation(lat=observatory.getLatitude().asDegrees()*units.deg,
620 lon=observatory.getLongitude().asDegrees()*units.deg,
621 height=observatory.getElevation()*units.m)
622 obstime = Time(visitInfo.getDate().get(system=DateTime.MJD),
623 location=loc, format="mjd")
624 coord = SkyCoord(
625 summary.ra*units.degree,
626 summary.dec*units.degree,
627 obstime=obstime,
628 location=loc,
629 )
630 with warnings.catch_warnings():
631 warnings.simplefilter("ignore")
632 altaz = coord.transform_to(AltAz)
634 summary.zenithDistance = float(90.0 - altaz.alt.degree)
636 def update_photo_calib_stats(self, summary, photo_calib):
637 """Compute all summary-statistic fields that depend on the photometric
638 calibration model.
640 Parameters
641 ----------
642 summary : `lsst.afw.image.ExposureSummaryStats`
643 Summary object to update in-place.
644 photo_calib : `lsst.afw.image.PhotoCalib` or `None`
645 Photometric calibration model. If `None`, all fields that depend
646 on the photometric calibration will be reset (generally to NaN).
647 """
648 if photo_calib is not None:
649 summary.zeroPoint = float(2.5*np.log10(photo_calib.getInstFluxAtZeroMagnitude()))
650 else:
651 summary.zeroPoint = float("nan")
653 def update_background_stats(self, summary, background):
654 """Compute summary-statistic fields that depend only on the
655 background model.
657 Parameters
658 ----------
659 summary : `lsst.afw.image.ExposureSummaryStats`
660 Summary object to update in-place.
661 background : `lsst.afw.math.BackgroundList` or `None`
662 Background model. If `None`, all fields that depend on the
663 background will be reset (generally to NaN).
665 Notes
666 -----
667 This does not include fields that depend on the background-subtracted
668 masked image; when the background changes, it should generally be
669 applied to the image and `update_masked_image_stats` should be called
670 as well.
671 """
672 if background is not None:
673 bgStats = (bg[0].getStatsImage().getImage().array
674 for bg in background)
675 summary.skyBg = float(sum(np.median(bg[np.isfinite(bg)]) for bg in bgStats))
676 else:
677 summary.skyBg = float("nan")
679 def update_masked_image_stats(self, summary, masked_image):
680 """Compute summary-statistic fields that depend on the masked image
681 itself.
683 Parameters
684 ----------
685 summary : `lsst.afw.image.ExposureSummaryStats`
686 Summary object to update in-place.
687 masked_image : `lsst.afw.image.MaskedImage` or `None`
688 Masked image. If `None`, all fields that depend
689 on the masked image will be reset (generally to NaN).
690 """
691 nan = float("nan")
692 if masked_image is None:
693 summary.skyNoise = nan
694 summary.meanVar = nan
695 return
696 statsCtrl = afwMath.StatisticsControl()
697 statsCtrl.setNumSigmaClip(self.config.sigmaClip)
698 statsCtrl.setNumIter(self.config.clipIter)
699 statsCtrl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.badMaskPlanes))
700 statsCtrl.setNanSafe(True)
702 statObj = afwMath.makeStatistics(masked_image, afwMath.STDEVCLIP, statsCtrl)
703 skyNoise, _ = statObj.getResult(afwMath.STDEVCLIP)
704 summary.skyNoise = skyNoise
706 statObj = afwMath.makeStatistics(masked_image.variance, masked_image.mask, afwMath.MEANCLIP,
707 statsCtrl)
708 meanVar, _ = statObj.getResult(afwMath.MEANCLIP)
709 summary.meanVar = meanVar
711 def update_effective_time_stats(self, summary, exposure):
712 """Compute effective exposure time statistics to estimate depth.
714 The effective exposure time is the equivalent shutter open
715 time that would be needed under nominal conditions to give the
716 same signal-to-noise for a point source as what is achieved by
717 the observation of interest. This metric combines measurements
718 of the point-spread function, the sky brightness, and the
719 transparency. It assumes that the observation is
720 sky-background dominated.
722 .. _teff_definitions:
724 The effective exposure time and its subcomponents are defined in [1]_.
726 References
727 ----------
729 .. [1] Neilsen, E.H., Bernstein, G., Gruendl, R., and Kent, S. (2016).
730 Limiting Magnitude, \tau, teff, and Image Quality in DES Year 1
731 https://www.osti.gov/biblio/1250877/
733 Parameters
734 ----------
735 summary : `lsst.afw.image.ExposureSummaryStats`
736 Summary object to update in-place.
737 exposure : `lsst.afw.image.ExposureF`
738 Exposure to grab band and exposure time metadata
740 """
741 nan = float("nan")
742 summary.effTime = nan
743 summary.effTimePsfSigmaScale = nan
744 summary.effTimeSkyBgScale = nan
745 summary.effTimeZeroPointScale = nan
747 exposureTime = exposure.getInfo().getVisitInfo().getExposureTime()
748 filterLabel = exposure.getFilter()
749 if (filterLabel is None) or (not filterLabel.hasBandLabel):
750 band = None
751 else:
752 band = filterLabel.bandLabel
754 if band is None:
755 self.log.warning("No band associated with exposure; effTime not calculated.")
756 return
758 # PSF component
759 if np.isnan(summary.psfSigma):
760 self.log.debug("PSF sigma is NaN")
761 f_eff = nan
762 elif band not in self.config.fiducialPsfSigma:
763 self.log.debug(f"Fiducial PSF value not found for {band}")
764 f_eff = nan
765 else:
766 fiducialPsfSigma = self.config.fiducialPsfSigma[band]
767 f_eff = (summary.psfSigma / fiducialPsfSigma)**-2
769 # Transparency component
770 # Note: Assumes that the zeropoint includes the exposure time
771 if np.isnan(summary.zeroPoint):
772 self.log.debug("Zero point is NaN")
773 c_eff = nan
774 elif band not in self.config.fiducialZeroPoint:
775 self.log.debug(f"Fiducial zero point value not found for {band}")
776 c_eff = nan
777 else:
778 fiducialZeroPoint = self.config.fiducialZeroPoint[band]
779 zeroPointDiff = fiducialZeroPoint - (summary.zeroPoint - 2.5*np.log10(exposureTime))
780 c_eff = min(10**(-2.0*(zeroPointDiff)/2.5), self.config.maxEffectiveTransparency)
782 # Sky brightness component (convert to cts/s)
783 if np.isnan(summary.skyBg):
784 self.log.debug("Sky background is NaN")
785 b_eff = nan
786 elif band not in self.config.fiducialSkyBackground:
787 self.log.debug(f"Fiducial sky background value not found for {band}")
788 b_eff = nan
789 else:
790 fiducialSkyBackground = self.config.fiducialSkyBackground[band]
791 b_eff = fiducialSkyBackground/(summary.skyBg/exposureTime)
793 # Old effective exposure time scale factor
794 # t_eff = f_eff * c_eff * b_eff
796 # New effective exposure time (seconds) from magnitude limit
797 if band not in self.config.fiducialMagLim:
798 self.log.debug(f"Fiducial magnitude limit not found for {band}")
799 effectiveTime = nan
800 elif band not in self.config.fiducialExpTime:
801 self.log.debug(f"Fiducial exposure time not found for {band}")
802 effectiveTime = nan
803 else:
804 effectiveTime = compute_effective_time(summary.magLim,
805 self.config.fiducialMagLim[band],
806 self.config.fiducialExpTime[band])
808 # Output quantities
809 summary.effTime = float(effectiveTime)
810 summary.effTimePsfSigmaScale = float(f_eff)
811 summary.effTimeSkyBgScale = float(b_eff)
812 summary.effTimeZeroPointScale = float(c_eff)
814 def update_magnitude_limit_stats(self, summary, exposure):
815 """Compute a summary statistic for the point-source magnitude limit at
816 a given signal-to-noise ratio using exposure-level metadata.
818 The magnitude limit is calculated at a given SNR from the PSF,
819 sky, zeropoint, and readnoise in accordance with SMTN-002 [1]_,
820 LSE-40 [2]_, and DMTN-296 [3]_.
822 References
823 ----------
825 .. [1] "Calculating LSST limiting magnitudes and SNR" (SMTN-002)
826 .. [2] "Photon Rates and SNR Calculations" (LSE-40)
827 .. [3] "Calculations of Image and Catalog Depth" (DMTN-296)
830 Parameters
831 ----------
832 summary : `lsst.afw.image.ExposureSummaryStats`
833 Summary object to update in-place.
834 exposure : `lsst.afw.image.ExposureF`
835 Exposure to grab band and exposure time metadata
836 """
837 if exposure.getDetector() is None:
838 summary.magLim = float("nan")
839 return
841 # Calculate the average readnoise [e-]
842 readNoiseList = list(ipIsr.getExposureReadNoises(exposure).values())
843 readNoise = np.nanmean(readNoiseList)
844 if np.isnan(readNoise):
845 readNoise = 0.0
846 self.log.warning("Read noise set to NaN! Setting readNoise to 0.0 to proceed.")
848 # Calculate the average gain [e-/ADU]
849 gainList = list(ipIsr.getExposureGains(exposure).values())
850 gain = np.nanmean(gainList)
851 if np.isnan(gain):
852 self.log.warning("Gain set to NaN! Setting magLim to NaN.")
853 summary.magLim = float("nan")
854 return
856 # Get the image units (default to "adu" if metadata key absent)
857 md = exposure.getMetadata()
858 if md.get("LSST ISR UNITS", "adu") == "electron":
859 gain = 1.0
861 # Convert readNoise to image units
862 readNoise /= gain
864 # Calculate the limiting magnitude.
865 # Note 1: Assumes that the image and readnoise have the same units
866 # Note 2: Assumes that the zeropoint includes the exposure time
867 magLim = compute_magnitude_limit(summary.psfArea, summary.skyBg,
868 summary.zeroPoint, readNoise,
869 gain, self.config.magLimSnr)
871 summary.magLim = float(magLim)
874def maximum_nearest_psf_distance(
875 image_mask,
876 psf_cat,
877 sampling=8,
878 bad_mask_bits=["BAD", "CR", "INTRP", "SAT", "SUSPECT", "NO_DATA", "EDGE"],
879):
880 """Compute the maximum distance of an unmasked pixel to its nearest PSF.
882 Parameters
883 ----------
884 image_mask : `lsst.afw.image.Mask`
885 The mask plane associated with the exposure.
886 psf_cat : `lsst.afw.table.SourceCatalog` or `astropy.table.Table`
887 Catalog containing only the stars used in the PSF modeling.
888 sampling : `int`
889 Sampling rate in each dimension to create the grid of points on which
890 to evaluate the distance to the nearest PSF star. The tradeoff is
891 between adequate sampling versus speed.
892 bad_mask_bits : `list` [`str`]
893 Mask bits required to be absent for a pixel to be considered
894 "unmasked".
896 Returns
897 -------
898 max_dist_to_nearest_psf : `float`
899 The maximum distance (in pixels) of an unmasked pixel to its nearest
900 PSF model star.
901 """
902 mask_arr = image_mask.array[::sampling, ::sampling]
903 bitmask = image_mask.getPlaneBitMask(bad_mask_bits)
904 good = ((mask_arr & bitmask) == 0)
906 x = np.arange(good.shape[1]) * sampling
907 y = np.arange(good.shape[0]) * sampling
908 xx, yy = np.meshgrid(x, y)
910 dist_to_nearest_psf = np.full(good.shape, np.inf)
911 for psf in psf_cat:
912 x_psf = psf["slot_Centroid_x"]
913 y_psf = psf["slot_Centroid_y"]
914 dist_to_nearest_psf = np.minimum(dist_to_nearest_psf, np.hypot(xx - x_psf, yy - y_psf))
915 unmasked_dists = dist_to_nearest_psf * good
916 max_dist_to_nearest_psf = np.max(unmasked_dists)
918 return max_dist_to_nearest_psf
921def compute_psf_image_deltas(
922 image_mask,
923 image_psf,
924 sampling=96,
925 ap_radius_pix=3.0,
926 bad_mask_bits=["BAD", "CR", "INTRP", "SAT", "SUSPECT", "NO_DATA", "EDGE"],
927):
928 """Compute the delta between the maximum and minimum model PSF trace radius
929 values evaluated on a grid of points lying in the unmasked region of the
930 image.
932 Parameters
933 ----------
934 image_mask : `lsst.afw.image.Mask`
935 The mask plane associated with the exposure.
936 image_psf : `lsst.afw.detection.Psf`
937 The PSF model associated with the exposure.
938 sampling : `int`, optional
939 Sampling rate in each dimension to create the grid of points at which
940 to evaluate ``image_psf``s trace radius value. The tradeoff is between
941 adequate sampling versus speed.
942 ap_radius_pix : `float`, optional
943 Radius in pixels of the aperture on which to measure the flux of the
944 PSF model.
945 bad_mask_bits : `list` [`str`], optional
946 Mask bits required to be absent for a pixel to be considered
947 "unmasked".
949 Returns
950 -------
951 psf_trace_radius_delta, psf_ap_flux_delta : `float`
952 The delta (in pixels) between the maximum and minimum model PSF trace
953 radius values and the PSF aperture fluxes (with aperture radius of
954 max(2, 3*psfSigma)) evaluated on the x,y-grid subsampled on the
955 unmasked detector pixels by a factor of ``sampling``. If both the
956 model PSF trace radius value and aperture flux value on the grid
957 evaluate to NaN, then NaNs are returned immediately.
958 """
959 psf_trace_radius_list = []
960 psf_ap_flux_list = []
961 mask_arr = image_mask.array[::sampling, ::sampling]
962 bitmask = image_mask.getPlaneBitMask(bad_mask_bits)
963 good = ((mask_arr & bitmask) == 0)
965 x = np.arange(good.shape[1]) * sampling
966 y = np.arange(good.shape[0]) * sampling
967 xx, yy = np.meshgrid(x, y)
969 for x_mesh, y_mesh, good_mesh in zip(xx, yy, good):
970 for x_point, y_point, is_good in zip(x_mesh, y_mesh, good_mesh):
971 if is_good:
972 position = geom.Point2D(x_point, y_point)
973 psf_trace_radius = image_psf.computeShape(position).getTraceRadius()
974 psf_ap_flux = image_psf.computeApertureFlux(ap_radius_pix, position)
975 if ~np.isfinite(psf_trace_radius) and ~np.isfinite(psf_ap_flux):
976 return float("nan"), float("nan")
977 psf_trace_radius_list.append(psf_trace_radius)
978 psf_ap_flux_list.append(psf_ap_flux)
980 psf_trace_radius_delta = np.max(psf_trace_radius_list) - np.min(psf_trace_radius_list)
981 if np.any(np.asarray(psf_ap_flux_list) < 0.0): # Consider any -ve flux entry as "bad".
982 psf_ap_flux_delta = float("nan")
983 else:
984 psf_ap_flux_delta = np.max(psf_ap_flux_list) - np.min(psf_ap_flux_list)
986 return psf_trace_radius_delta, psf_ap_flux_delta
989def compute_ap_corr_sigma_scaled_delta(
990 image_mask,
991 image_ap_corr_field,
992 psfSigma,
993 sampling=96,
994 bad_mask_bits=["BAD", "CR", "INTRP", "SAT", "SUSPECT", "NO_DATA", "EDGE"],
995):
996 """Compute the delta between the maximum and minimum aperture correction
997 values scaled (divided) by ``psfSigma`` for the given field representation,
998 ``image_ap_corr_field`` evaluated on a grid of points lying in the
999 unmasked region of the image.
1001 Parameters
1002 ----------
1003 image_mask : `lsst.afw.image.Mask`
1004 The mask plane associated with the exposure.
1005 image_ap_corr_field : `lsst.afw.math.ChebyshevBoundedField`
1006 The ChebyshevBoundedField representation of the aperture correction
1007 of interest for the exposure.
1008 psfSigma : `float`
1009 The PSF model second-moments determinant radius (center of chip)
1010 in pixels.
1011 sampling : `int`, optional
1012 Sampling rate in each dimension to create the grid of points at which
1013 to evaluate ``image_psf``s trace radius value. The tradeoff is between
1014 adequate sampling versus speed.
1015 bad_mask_bits : `list` [`str`], optional
1016 Mask bits required to be absent for a pixel to be considered
1017 "unmasked".
1019 Returns
1020 -------
1021 ap_corr_sigma_scaled_delta : `float`
1022 The delta between the maximum and minimum of the (multiplicative)
1023 aperture correction values scaled (divided) by ``psfSigma`` evaluated
1024 on the x,y-grid subsampled on the unmasked detector pixels by a factor
1025 of ``sampling``. If the aperture correction evaluates to NaN on any
1026 of the grid points, this is set to NaN.
1027 """
1028 mask_arr = image_mask.array[::sampling, ::sampling]
1029 bitmask = image_mask.getPlaneBitMask(bad_mask_bits)
1030 good = ((mask_arr & bitmask) == 0)
1032 x = np.arange(good.shape[1], dtype=np.float64) * sampling
1033 y = np.arange(good.shape[0], dtype=np.float64) * sampling
1034 xx, yy = np.meshgrid(x, y)
1036 ap_corr = image_ap_corr_field.evaluate(xx.ravel(), yy.ravel()).reshape(xx.shape)
1037 ap_corr_good = ap_corr[good]
1038 if ~np.isfinite(ap_corr_good).all():
1039 ap_corr_sigma_scaled_delta = float("nan")
1040 else:
1041 ap_corr_sigma_scaled_delta = (np.max(ap_corr_good) - np.min(ap_corr_good))/psfSigma
1043 return ap_corr_sigma_scaled_delta
1046def compute_magnitude_limit(
1047 psfArea,
1048 skyBg,
1049 zeroPoint,
1050 readNoise,
1051 gain,
1052 snr
1053):
1054 """Compute the expected point-source magnitude limit at a given
1055 signal-to-noise ratio given the exposure-level metadata. Based on
1056 the signal-to-noise formula provided in SMTN-002 (see LSE-40 for
1057 more details on the calculation).
1059 SNR = C / sqrt( C/g + (B/g + sigma_inst**2) * neff )
1061 where C is the counts from the source, B is counts from the (sky)
1062 background, sigma_inst is the instrumental (read) noise, neff is
1063 the effective size of the PSF, and g is the gain in e-/ADU. Note
1064 that input values of ``skyBg``, ``zeroPoint``, and ``readNoise``
1065 should all consistently be in electrons or ADU.
1067 Parameters
1068 ----------
1069 psfArea : `float`
1070 The effective area of the PSF [pix].
1071 skyBg : `float`
1072 The sky background counts for the exposure [ADU or e-].
1073 zeroPoint : `float`
1074 The zeropoint (includes exposure time) [ADU or e-].
1075 readNoise : `float`
1076 The instrumental read noise for the exposure [ADU or e-].
1077 gain : `float`
1078 The instrumental gain for the exposure [e-/ADU]. The gain should
1079 be 1.0 if the skyBg, zeroPoint, and readNoise are in e-.
1080 snr : `float`
1081 Signal-to-noise ratio at which magnitude limit is calculated.
1083 Returns
1084 -------
1085 magnitude_limit : `float`
1086 The expected magnitude limit at the given signal to noise.
1087 """
1088 # Effective number of pixels within the PSF calculated directly
1089 # from the PSF model.
1090 neff = psfArea
1092 # Instrumental noise (read noise only)
1093 sigma_inst = readNoise
1095 # Total background counts derived from Eq. 45 of LSE-40
1096 background = (skyBg/gain + sigma_inst**2) * neff
1097 # Source counts to achieve the desired SNR (from quadratic formula)
1098 # Note typo in gain exponent in Eq. 45 of LSE-40 (DM-53234)
1099 source = (snr**2)/(2*gain) + np.sqrt((snr**4)/(4*gain**2) + snr**2 * background)
1101 # Convert to a magnitude using the zeropoint.
1102 # Note: Zeropoint currently includes exposure time
1103 magnitude_limit = -2.5*np.log10(source) + zeroPoint
1105 return magnitude_limit
1108def psf_sigma_to_psf_area(psfSigma, pixelScale):
1109 """Convert psfSigma [pixels] to an approximate psfArea [pixel^2] as defined in LSE-40.
1110 This is the same implementation followed by SMTN-002 to estimate SNR=5 magnitude limits.
1112 Parameters
1113 ----------
1114 psfSigma : `float`
1115 The PSF sigma value [pix].
1116 pixelScale : `float`
1117 The pixel scale [arcsec/pix].
1119 Returns
1120 -------
1121 psf_area : `float`
1122 Approximation of the PSF area [pix^2]
1123 """
1124 # Follow SMTN-002 to convert to geometric and effective FWHM
1125 fwhm_geom = psfSigma * 2.355 * pixelScale
1126 fwhm_eff = (fwhm_geom - 0.052) / 0.822
1127 # Area prefactor comes from LSE-40
1128 psf_area = 2.266 * (fwhm_eff / pixelScale)**2
1129 return psf_area
1132def compute_effective_time(
1133 magLim,
1134 fiducialMagLim,
1135 fiducialExpTime
1136):
1137 """Compute the effective exposure time from m5 following the prescription described in SMTN-296.
1139 teff = 10**(0.8 * (magLim - fiducialMagLim) ) * fiducialExpTime
1141 where `magLim` is the magnitude limit, `fiducialMagLim` is the fiducial magnitude limit calculated from
1142 the fiducial values of the ``psfArea``, ``skyBg``, ``zeroPoint``, and ``readNoise``, and
1143 `fiducialExpTime` is the fiducial exposure time (s).
1145 Parameters
1146 ----------
1147 magLim : `float`
1148 The measured magnitude limit [mag].
1149 fiducialMagLim : `float`
1150 The fiducial magnitude limit [mag].
1151 fiducialExpTime : `float`
1152 The fiducial exposure time [s].
1154 Returns
1155 -------
1156 effectiveTime : `float`
1157 The effective exposure time.
1158 """
1159 effectiveTime = 10**(0.8 * (magLim - fiducialMagLim)) * fiducialExpTime
1160 return effectiveTime