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