Coverage for python/lsst/pipe/tasks/computeExposureSummaryStats.py: 72%
528 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-23 09: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.meas.algorithms.computeExPsf import ComputeExPsfTask
39from lsst.utils.timer import timeMethod
40import lsst.ip.isr as ipIsr
43class ComputeExposureSummaryStatsConfig(pexConfig.Config):
44 """Config for ComputeExposureSummaryTask"""
45 sigmaClip = pexConfig.Field(
46 dtype=float,
47 doc="Sigma for outlier rejection for sky noise.",
48 default=3.0,
49 )
50 clipIter = pexConfig.Field(
51 dtype=int,
52 doc="Number of iterations of outlier rejection for sky noise.",
53 default=2,
54 )
55 badMaskPlanes = pexConfig.ListField(
56 dtype=str,
57 doc="Mask planes that, if set, the associated pixel should not be included sky noise calculation.",
58 default=("NO_DATA", "SUSPECT"),
59 )
60 starSelection = pexConfig.Field(
61 doc="Field to select full list of sources used for PSF modeling.",
62 dtype=str,
63 default="calib_psf_used",
64 )
65 starSelector = pexConfig.ConfigurableField(
66 target=ScienceSourceSelectorTask,
67 doc="Selection of sources to compute PSF star statistics.",
68 )
69 starShape = pexConfig.Field(
70 doc="Base name of columns to use for the source shape in the PSF statistics computation.",
71 dtype=str,
72 default="slot_Shape"
73 )
74 psfShape = pexConfig.Field(
75 doc="Base name of columns to use for the PSF shape in the PSF statistics computation.",
76 dtype=str,
77 default="slot_PsfShape"
78 )
79 starHigherOrderMomentBase = pexConfig.Field(
80 doc="Base name of the columns from which to obtain the higher-order moments.",
81 dtype=str,
82 default="ext_shapeHSM_HigherOrderMomentsSource",
83 )
84 psfSampling = pexConfig.Field(
85 dtype=int,
86 doc="Sampling rate in pixels in each dimension for the maxDistToNearestPsf metric "
87 "caclulation grid (the tradeoff is between adequate sampling versus speed).",
88 default=8,
89 )
90 psfGridSampling = pexConfig.Field(
91 dtype=int,
92 doc="Sampling rate in pixels in each dimension for PSF model robustness metric "
93 "caclulations grid (the tradeoff is between adequate sampling versus speed).",
94 default=96,
95 )
96 minPsfApRadiusPix = pexConfig.Field(
97 dtype=float,
98 doc="Minimum radius in pixels of the aperture within which to measure the flux of "
99 "the PSF model for the psfApFluxDelta metric calculation (the radius is computed as "
100 "max(``minPsfApRadius``, 3*psfSigma)).",
101 default=2.0,
102 )
103 psfApCorrFieldName = pexConfig.Field(
104 doc="Name of the flux column associated with the aperture correction of the PSF model "
105 "to use for the psfApCorrSigmaScaledDelta metric calculation.",
106 dtype=str,
107 default="base_PsfFlux_instFlux"
108 )
109 psfBadMaskPlanes = pexConfig.ListField(
110 dtype=str,
111 doc="Mask planes that, if set, the associated pixel should not be included in the PSF model "
112 "robutsness metric calculations (namely, maxDistToNearestPsf and psfTraceRadiusDelta).",
113 default=("BAD", "CR", "EDGE", "INTRP", "NO_DATA", "SAT", "SUSPECT"),
114 )
115 fiducialSkyBackground = pexConfig.DictField(
116 keytype=str,
117 itemtype=float,
118 doc="Fiducial sky background level (ADU/s) assumed when calculating effective exposure time. "
119 "Keyed by band.",
120 default={"u": 1.0, "g": 1.0, "r": 1.0, "i": 1.0, "z": 1.0, "y": 1.0},
121 )
122 fiducialPsfSigma = pexConfig.DictField(
123 keytype=str,
124 itemtype=float,
125 doc="Fiducial PSF sigma (pixels) assumed when calculating effective exposure time. "
126 "Keyed by band.",
127 default={"u": 1.0, "g": 1.0, "r": 1.0, "i": 1.0, "z": 1.0, "y": 1.0},
128 )
129 fiducialZeroPoint = pexConfig.DictField(
130 keytype=str,
131 itemtype=float,
132 doc="Fiducial zero point assumed when calculating effective exposure time. "
133 "Keyed by band.",
134 default={"u": 25.0, "g": 25.0, "r": 25.0, "i": 25.0, "z": 25.0, "y": 25.0},
135 )
136 fiducialReadNoise = pexConfig.DictField(
137 keytype=str,
138 itemtype=float,
139 doc="Fiducial readnoise (electrons) assumed when calculating effective exposure time. "
140 "Keyed by band.",
141 default={"u": 9.0, "g": 9.0, "r": 9.0, "i": 9.0, "z": 9.0, "y": 9.0},
142 )
143 fiducialExpTime = pexConfig.DictField(
144 keytype=str,
145 itemtype=float,
146 doc="Fiducial exposure time (seconds). "
147 "Keyed by band.",
148 default={"u": 30.0, "g": 30.0, "r": 30.0, "i": 30.0, "z": 30.0, "y": 30.0},
149 )
150 fiducialMagLim = pexConfig.DictField(
151 keytype=str,
152 itemtype=float,
153 doc="Fiducial magnitude limit depth at SNR=5. "
154 "Keyed by band.",
155 default={"u": 25.0, "g": 25.0, "r": 25.0, "i": 25.0, "z": 25.0, "y": 25.0},
156 )
157 maxEffectiveTransparency = pexConfig.Field(
158 dtype=float,
159 doc="Maximum value allowed for effective transparency scale factor (often inf or 1.0).",
160 default=float("inf")
161 )
162 magLimSnr = pexConfig.Field(
163 dtype=float,
164 doc="Signal-to-noise ratio for computing the magnitude limit depth.",
165 default=5.0
166 )
167 psfTE1 = pexConfig.ConfigurableField(
168 target=ComputeExPsfTask,
169 doc="Use treecorr for computing scalar value of TE1.",
170 )
171 psfTE2 = pexConfig.ConfigurableField(
172 target=ComputeExPsfTask,
173 doc="Use treecorr for computing scalar value of TE2.",
174 )
175 psfTE3 = pexConfig.ConfigurableField(
176 target=ComputeExPsfTask,
177 doc="Use treecorr for computing scalar value of TE3.",
178 )
179 psfTE4 = pexConfig.ConfigurableField(
180 target=ComputeExPsfTask,
181 doc="Use treecorr for computing scalar value of TE4.",
182 )
184 def setDefaults(self):
185 super().setDefaults()
187 self.starSelector.setDefaults()
188 self.starSelector.doFlags = True
189 self.starSelector.doSignalToNoise = True
190 self.starSelector.doUnresolved = False
191 self.starSelector.doIsolated = False
192 self.starSelector.doRequireFiniteRaDec = False
193 self.starSelector.doRequirePrimary = False
195 self.starSelector.signalToNoise.minimum = 50.0
196 self.starSelector.signalToNoise.maximum = 1000.0
198 self.starSelector.flags.bad = ["slot_Shape_flag", "slot_PsfFlux_flag"]
199 # Select stars used for PSF modeling.
200 self.starSelector.flags.good = ["calib_psf_used"]
202 self.starSelector.signalToNoise.fluxField = "slot_PsfFlux_instFlux"
203 self.starSelector.signalToNoise.errField = "slot_PsfFlux_instFluxErr"
205 min_theta = [1e-6, 5.0, 1e-6, 5.0]
206 max_theta = [1.0, 100.0, 5.0, 20.0]
207 psfTEx = [
208 self.psfTE1,
209 self.psfTE2,
210 self.psfTE3,
211 self.psfTE4,
212 ]
214 for tex, mint, maxt in zip(psfTEx, min_theta, max_theta):
215 tex.setDefaults()
216 tex.treecorr.min_sep = mint / 60.0
217 tex.treecorr.max_sep = maxt / 60.0
218 tex.treecorr.nbins = 1
219 tex.treecorr.bin_type = "Linear"
220 tex.treecorr.sep_units = "degree"
222 def fiducialMagnitudeLimit(self, band, pixelScale, gain):
223 """Compute the fiducial point-source magnitude limit based on config values.
224 This follows the conventions laid out in SMTN-002, LSE-40, and DMTN-296.
226 Parameters
227 ----------
228 band : `str`
229 The band of interest
230 pixelScale : `float`
231 The pixel scale [arcsec/pix]
232 gain : `float`
233 The instrumental gain for the exposure [e-/ADU]. The gain should
234 be 1.0 if the image units are [e-].
236 Returns
237 -------
238 magnitude_limit : `float`
239 The fiducial magnitude limit calculated from fiducial values.
240 """
241 nan = float("nan")
243 # Fiducials from config
244 fiducialPsfSigma = self.fiducialPsfSigma.get(band, nan)
245 fiducialSkyBackground = self.fiducialSkyBackground.get(band, nan)
246 fiducialZeroPoint = self.fiducialZeroPoint.get(band, nan)
247 fiducialReadNoise = self.fiducialReadNoise.get(band, nan)
248 fiducialExpTime = self.fiducialExpTime.get(band, nan)
249 magLimSnr = self.magLimSnr
251 # Derived fiducial quantities
252 fiducialPsfArea = psf_sigma_to_psf_area(fiducialPsfSigma, pixelScale)
253 fiducialZeroPoint += 2.5*np.log10(fiducialExpTime)
254 fiducialSkyBg = fiducialSkyBackground * fiducialExpTime
255 fiducialReadNoise /= gain
257 # Calculate the fiducial magnitude limit
258 fiducialMagLim = compute_magnitude_limit(fiducialPsfArea,
259 fiducialSkyBg,
260 fiducialZeroPoint,
261 fiducialReadNoise,
262 gain,
263 magLimSnr)
265 return fiducialMagLim
268class ComputeExposureSummaryStatsTask(pipeBase.Task):
269 """Task to compute exposure summary statistics.
271 This task computes various quantities suitable for DPDD and other
272 downstream processing at the detector centers, including:
273 - expTime
274 - psfSigma
275 - psfArea
276 - psfIxx
277 - psfIyy
278 - psfIxy
279 - ra
280 - dec
281 - pixelScale (arcsec/pixel)
282 - zenithDistance
283 - zeroPoint
284 - skyBg
285 - skyNoise
286 - meanVar
287 - raCorners
288 - decCorners
289 - astromOffsetMean
290 - astromOffsetStd
292 These additional quantities are computed from the stars in the detector:
293 - psfStarDeltaE1Median
294 - psfStarDeltaE2Median
295 - psfStarDeltaE1Scatter
296 - psfStarDeltaE2Scatter
297 - psfStarDeltaSizeMedian
298 - psfStarDeltaSizeScatter
299 - psfStarScaledDeltaSizeScatter
301 These quantities are computed based on the PSF model and image mask
302 to assess the robustness of the PSF model across a given detector
303 (against, e.g., extrapolation instability):
304 - maxDistToNearestPsf
305 - psfTraceRadiusDelta
306 - psfApFluxDelta
308 These quantities are computed as part of:
309 https://rubinobs.atlassian.net/browse/DM-40780
311 - psfTE1e1
312 - psfTE1e2
313 - psfTE1ex
314 - psfTE2e1
315 - psfTE2e2
316 - psfTE2ex
317 - psfTE3e1
318 - psfTE3e2
319 - psfTE3ex
320 - psfTE4e1
321 - psfTE4e2
322 - psfTE4ex
324 This quantity is computed based on the aperture correction map, the
325 psfSigma, and the image mask to assess the robustness of the aperture
326 corrections across a given detector:
327 - psfApCorrSigmaScaledDelta
329 These quantities are computed based on the shape measurements of the
330 sources used in the PSF model to assess the image quality (i.e. as a
331 measure of the deviation from a circular shape):
332 - starEMedian
333 - starUnNormalizedEMedian
335 These quantities are computed to assess depth:
336 - effTime
337 - effTimePsfSigmaScale
338 - effTimeSkyBgScale
339 - effTimeZeroPointScale
340 - magLim
341 """
342 ConfigClass = ComputeExposureSummaryStatsConfig
343 _DefaultName = "computeExposureSummaryStats"
345 def __init__(self, **kwargs):
346 super().__init__(**kwargs)
348 self.makeSubtask("starSelector")
349 self.makeSubtask("psfTE1")
350 self.makeSubtask("psfTE2")
351 self.makeSubtask("psfTE3")
352 self.makeSubtask("psfTE4")
353 self._isTEXComputationDone = False
355 @timeMethod
356 def run(self, exposure, sources, background, summary=None):
357 """Measure exposure statistics from the exposure, sources, and
358 background.
360 Parameters
361 ----------
362 exposure : `lsst.afw.image.ExposureF`
363 sources : `lsst.afw.table.SourceCatalog`
364 background : `lsst.afw.math.BackgroundList`
365 summary : `lsst.afw.image.ExposureSummary`, optional
366 Summary object to be appended to if not `None`. If `None` a new
367 summary object will be created.
369 Returns
370 -------
371 summary : `lsst.afw.image.ExposureSummary`
372 """
373 self.log.info("Measuring exposure statistics")
375 if summary is None:
376 summary = afwImage.ExposureSummaryStats()
378 # Set exposure time.
379 exposureTime = exposure.getInfo().getVisitInfo().getExposureTime()
380 summary.expTime = exposureTime
382 bbox = exposure.getBBox()
384 psf = exposure.getPsf()
385 self.update_psf_stats(
386 summary, psf, bbox, sources, image_mask=exposure.mask, image_ap_corr_map=exposure.apCorrMap
387 )
389 wcs = exposure.getWcs()
390 visitInfo = exposure.getInfo().getVisitInfo()
391 self.update_wcs_stats(summary, wcs, bbox, visitInfo)
393 photoCalib = exposure.getPhotoCalib()
394 self.update_photo_calib_stats(summary, photoCalib)
396 self.update_background_stats(summary, background)
398 self.update_masked_image_stats(summary, exposure.getMaskedImage())
400 self.update_magnitude_limit_stats(summary, exposure)
402 self.update_effective_time_stats(summary, exposure)
404 md = exposure.getMetadata()
405 if "PSF_ADAPTIVE_THRESHOLD_VALUE" in md:
406 summary.psfAdaptiveThresholdValue = md["PSF_ADAPTIVE_THRESHOLD_VALUE"]
407 if "PSF_ADAPTIVE_INCLUDE_THRESHOLD_MULTIPLIER" in md:
408 summary.psfAdaptiveIncludeThresholdMultiplier = md["PSF_ADAPTIVE_INCLUDE_THRESHOLD_MULTIPLIER"]
409 if "REF_CAT_SOURCE_DENSITY" in md:
410 summary.refCatSourceDensity = md["REF_CAT_SOURCE_DENSITY"]
411 if "SFM_ASTROM_OFFSET_MEAN" in md:
412 summary.astromOffsetMean = md["SFM_ASTROM_OFFSET_MEAN"]
413 summary.astromOffsetStd = md["SFM_ASTROM_OFFSET_STD"]
415 return summary
417 def update_psf_stats(
418 self,
419 summary,
420 psf,
421 bbox,
422 sources=None,
423 image_mask=None,
424 image_ap_corr_map=None,
425 sources_is_astropy=False,
426 ):
427 """Compute all summary-statistic fields that depend on the PSF model.
429 Parameters
430 ----------
431 summary : `lsst.afw.image.ExposureSummaryStats`
432 Summary object to update in-place.
433 psf : `lsst.afw.detection.Psf` or `None`
434 Point spread function model. If `None`, all fields that depend on
435 the PSF will be reset (generally to NaN).
436 bbox : `lsst.geom.Box2I`
437 Bounding box of the image for which summary stats are being
438 computed.
439 sources : `lsst.afw.table.SourceCatalog` or `astropy.table.Table`
440 Catalog for quantities that are computed from source table columns.
441 If `None`, these quantities will be reset (generally to NaN).
442 The type of this table must correspond to the
443 ``sources_is_astropy`` argument.
444 image_mask : `lsst.afw.image.Mask`, optional
445 Mask image that may be used to compute distance-to-nearest-star
446 metrics.
447 sources_is_astropy : `bool`, optional
448 Whether ``sources`` is an `astropy.table.Table` instance instead
449 of an `lsst.afw.table.Catalog` instance. Default is `False` (the
450 latter).
451 """
452 nan = float("nan")
453 summary.psfSigma = nan
454 summary.psfIxx = nan
455 summary.psfIyy = nan
456 summary.psfIxy = nan
457 summary.psfArea = nan
458 summary.nPsfStar = 0
459 summary.psfStarDeltaE1Median = nan
460 summary.psfStarDeltaE2Median = nan
461 summary.psfStarDeltaE1Scatter = nan
462 summary.psfStarDeltaE2Scatter = nan
463 summary.psfStarDeltaSizeMedian = nan
464 summary.psfStarDeltaSizeScatter = nan
465 summary.psfStarScaledDeltaSizeScatter = nan
466 summary.maxDistToNearestPsf = nan
467 summary.psfTraceRadiusDelta = nan
468 summary.psfApFluxDelta = nan
469 summary.psfApCorrSigmaScaledDelta = nan
470 summary.starEMedian = nan
471 summary.starUnNormalizedEMedian = nan
472 summary.starComa1Median = nan
473 summary.starComa2Median = nan
474 summary.starTrefoil1Median = nan
475 summary.starTrefoil2Median = nan
476 summary.starKurtosisMedian = nan
477 summary.starE41Median = nan
478 summary.starE42Median = nan
480 if psf is None: 480 ↛ 481line 480 didn't jump to line 481 because the condition on line 480 was never true
481 return
482 shape = psf.computeShape(bbox.getCenter())
483 summary.psfSigma = shape.getDeterminantRadius()
484 summary.psfIxx = shape.getIxx()
485 summary.psfIyy = shape.getIyy()
486 summary.psfIxy = shape.getIxy()
487 im = psf.computeKernelImage(bbox.getCenter())
488 # The calculation of effective psf area is taken from
489 # ls.st/srd Equation 1.
490 summary.psfArea = float(np.sum(im.array)**2./np.sum(im.array**2.))
492 if image_mask is not None: 492 ↛ 520line 492 didn't jump to line 520 because the condition on line 492 was always true
493 psfApRadius = max(self.config.minPsfApRadiusPix, 3.0*summary.psfSigma)
494 self.log.debug("Using radius of %.3f (pixels) for psfApFluxDelta metric", psfApRadius)
495 psfTraceRadiusDelta, psfApFluxDelta = compute_psf_image_deltas(
496 image_mask,
497 psf,
498 sampling=self.config.psfGridSampling,
499 ap_radius_pix=psfApRadius,
500 bad_mask_bits=self.config.psfBadMaskPlanes
501 )
502 summary.psfTraceRadiusDelta = float(psfTraceRadiusDelta)
503 summary.psfApFluxDelta = float(psfApFluxDelta)
504 if image_ap_corr_map is not None:
505 if self.config.psfApCorrFieldName not in image_ap_corr_map.keys(): 505 ↛ 506line 505 didn't jump to line 506 because the condition on line 505 was never true
506 self.log.warning(f"{self.config.psfApCorrFieldName} not found in "
507 "image_ap_corr_map. Setting psfApCorrSigmaScaledDelta to NaN.")
508 psfApCorrSigmaScaledDelta = nan
509 else:
510 image_ap_corr_field = image_ap_corr_map[self.config.psfApCorrFieldName]
511 psfApCorrSigmaScaledDelta = compute_ap_corr_sigma_scaled_delta(
512 image_mask,
513 image_ap_corr_field,
514 summary.psfSigma,
515 sampling=self.config.psfGridSampling,
516 bad_mask_bits=self.config.psfBadMaskPlanes,
517 )
518 summary.psfApCorrSigmaScaledDelta = float(psfApCorrSigmaScaledDelta)
520 if sources is None:
521 # No sources are available (as in some tests and rare cases where
522 # the selection criteria in finalizeCharacterization lead to no
523 # good sources).
524 return
526 # Count the total number of psf stars used (prior to stats selection).
527 nPsfStar = sources[self.config.starSelection].sum()
528 summary.nPsfStar = int(nPsfStar)
530 psf_mask = self.starSelector.run(sources).selected
531 nPsfStarsUsedInStats = psf_mask.sum()
533 if nPsfStarsUsedInStats == 0: 533 ↛ 536line 533 didn't jump to line 536 because the condition on line 533 was never true
534 # No stars to measure statistics, so we must return the defaults
535 # of 0 stars and NaN values.
536 return
538 if sources_is_astropy: 538 ↛ 539line 538 didn't jump to line 539 because the condition on line 538 was never true
539 psf_cat = sources[psf_mask]
540 else:
541 psf_cat = sources[psf_mask].copy(deep=True)
543 starXX = psf_cat[self.config.starShape + "_xx"]
544 starYY = psf_cat[self.config.starShape + "_yy"]
545 starXY = psf_cat[self.config.starShape + "_xy"]
546 psfXX = psf_cat[self.config.psfShape + "_xx"]
547 psfYY = psf_cat[self.config.psfShape + "_yy"]
548 psfXY = psf_cat[self.config.psfShape + "_xy"]
550 # Use the trace radius for the star size.
551 starSize = np.sqrt(starXX/2. + starYY/2.)
553 starE1 = (starXX - starYY)/(starXX + starYY)
554 starE2 = 2*starXY/(starXX + starYY)
555 starSizeMedian = np.median(starSize)
557 starE = np.sqrt(starE1**2.0 + starE2**2.0)
558 starUnNormalizedE = np.sqrt((starXX - starYY)**2.0 + (2.0*starXY)**2.0)
559 starEMedian = np.median(starE)
560 starUnNormalizedEMedian = np.median(starUnNormalizedE)
561 summary.starEMedian = float(starEMedian)
562 summary.starUnNormalizedEMedian = float(starUnNormalizedEMedian)
564 # Use the trace radius for the psf size.
565 psfSize = np.sqrt(psfXX/2. + psfYY/2.)
566 psfE1 = (psfXX - psfYY)/(psfXX + psfYY)
567 psfE2 = 2*psfXY/(psfXX + psfYY)
569 e1Residuals = starE1 - psfE1
570 e2Residuals = starE2 - psfE2
572 psfStarDeltaE1Median = np.median(e1Residuals)
573 psfStarDeltaE1Scatter = sigmaMad(e1Residuals, scale='normal')
574 psfStarDeltaE2Median = np.median(e2Residuals)
575 psfStarDeltaE2Scatter = sigmaMad(e2Residuals, scale='normal')
577 psfStarDeltaSizeMedian = np.median(starSize - psfSize)
578 psfStarDeltaSizeScatter = sigmaMad(starSize - psfSize, scale="normal")
579 psfStarScaledDeltaSizeScatter = psfStarDeltaSizeScatter/starSizeMedian
581 summary.psfStarDeltaE1Median = float(psfStarDeltaE1Median)
582 summary.psfStarDeltaE2Median = float(psfStarDeltaE2Median)
583 summary.psfStarDeltaE1Scatter = float(psfStarDeltaE1Scatter)
584 summary.psfStarDeltaE2Scatter = float(psfStarDeltaE2Scatter)
585 summary.psfStarDeltaSizeMedian = float(psfStarDeltaSizeMedian)
586 summary.psfStarDeltaSizeScatter = float(psfStarDeltaSizeScatter)
587 summary.psfStarScaledDeltaSizeScatter = float(psfStarScaledDeltaSizeScatter)
589 # Higher-order moments and derived metrics.
590 if sources_is_astropy: 590 ↛ 591line 590 didn't jump to line 591 because the condition on line 590 was never true
591 columnNames = psf_cat.colnames
592 else:
593 columnNames = psf_cat.schema.getNames()
594 if self.config.starHigherOrderMomentBase + "_03" in columnNames: 594 ↛ 595line 594 didn't jump to line 595 because the condition on line 594 was never true
595 starM03 = psf_cat[self.config.starHigherOrderMomentBase + "_03"]
596 starM12 = psf_cat[self.config.starHigherOrderMomentBase + "_12"]
597 starM21 = psf_cat[self.config.starHigherOrderMomentBase + "_21"]
598 starM30 = psf_cat[self.config.starHigherOrderMomentBase + "_30"]
599 starM04 = psf_cat[self.config.starHigherOrderMomentBase + "_04"]
600 starM13 = psf_cat[self.config.starHigherOrderMomentBase + "_13"]
601 starM22 = psf_cat[self.config.starHigherOrderMomentBase + "_22"]
602 starM31 = psf_cat[self.config.starHigherOrderMomentBase + "_31"]
603 starM40 = psf_cat[self.config.starHigherOrderMomentBase + "_40"]
605 starComa1Median = np.nanmedian(starM30 + starM12)
606 starComa2Median = np.nanmedian(starM21 + starM03)
607 starTrefoil1Median = np.nanmedian(starM30 - 3.0*starM12)
608 starTrefoil2Median = np.nanmedian(3.0*starM21 - starM03)
609 starKurtosisMedian = np.nanmedian(starM40 + 2.0*starM22 + starM04)
610 starE41Median = np.nanmedian(starM40 - starM04)
611 starE42Median = np.nanmedian(2.0*(starM31 + starM13))
613 summary.starComa1Median = float(starComa1Median)
614 summary.starComa2Median = float(starComa2Median)
615 summary.starTrefoil1Median = float(starTrefoil1Median)
616 summary.starTrefoil2Median = float(starTrefoil2Median)
617 summary.starKurtosisMedian = float(starKurtosisMedian)
618 summary.starE41Median = float(starE41Median)
619 summary.starE42Median = float(starE42Median)
620 else:
621 self.log.warning("Higher-order moments with base column name %s not found in source "
622 "catalog. Setting the derived metrics (i.e. coma1[2], trefoil1[2], "
623 "kurtosis, e4_1, and e4_2) to nan.", self.config.starHigherOrderMomentBase)
625 if image_mask is not None: 625 ↛ exitline 625 didn't return from function 'update_psf_stats' because the condition on line 625 was always true
626 maxDistToNearestPsf = maximum_nearest_psf_distance(
627 image_mask,
628 psf_cat,
629 sampling=self.config.psfSampling,
630 bad_mask_bits=self.config.psfBadMaskPlanes
631 )
632 summary.maxDistToNearestPsf = float(maxDistToNearestPsf)
634 def comp_psf_TEX_visit_level(self, summary, sources, sources_is_astropy=False):
635 """Compute all summary-statistic fields at visit level for TEx metric of PSF.
637 Parameters
638 ----------
639 summary : `lsst.afw.image.ExposureSummaryStats`
640 Summary object to update in-place.
641 sources : `lsst.afw.table.SourceCatalog` or `astropy.table.Table`
642 Catalog for quantities that are computed from source table columns.
643 If `None`, these quantities will be reset (generally to NaN).
644 The type of this table must correspond to the
645 ``sources_is_astropy`` argument.
646 sources_is_astropy : `bool`, optional
647 Whether ``sources`` is an `astropy.table.Table` instance instead
648 of an `lsst.afw.table.Catalog` instance.
649 """
651 if self._isTEXComputationDone:
653 summary.psfTE1e1 = self.psfTE1e1
654 summary.psfTE1e2 = self.psfTE1e2
655 summary.psfTE1ex = self.psfTE1ex
656 summary.psfTE2e1 = self.psfTE2e1
657 summary.psfTE2e2 = self.psfTE2e2
658 summary.psfTE2ex = self.psfTE2ex
659 summary.psfTE3e1 = self.psfTE3e1
660 summary.psfTE3e2 = self.psfTE3e2
661 summary.psfTE3ex = self.psfTE3ex
662 summary.psfTE4e1 = self.psfTE4e1
663 summary.psfTE4e2 = self.psfTE4e2
664 summary.psfTE4ex = self.psfTE4ex
666 else:
668 self._isTEXComputationDone = True
670 nan = float("nan")
671 summary.psfTE1e1, self.psfTE1e1 = nan, nan
672 summary.psfTE1e2, self.psfTE1e2 = nan, nan
673 summary.psfTE1ex, self.psfTE1ex = nan, nan
674 summary.psfTE2e1, self.psfTE2e1 = nan, nan
675 summary.psfTE2e2, self.psfTE2e2 = nan, nan
676 summary.psfTE2ex, self.psfTE2ex = nan, nan
677 summary.psfTE3e1, self.psfTE3e1 = nan, nan
678 summary.psfTE3e2, self.psfTE3e2 = nan, nan
679 summary.psfTE3ex, self.psfTE3ex = nan, nan
680 summary.psfTE4e1, self.psfTE4e1 = nan, nan
681 summary.psfTE4e2, self.psfTE4e2 = nan, nan
682 summary.psfTE4ex, self.psfTE4ex = nan, nan
684 psf_mask = self.starSelector.run(sources).selected
686 nPsfStarsUsedInStats = psf_mask.sum()
688 if nPsfStarsUsedInStats == 0:
689 # No stars to measure statistics, so we must return the defaults
690 # of 0 stars and NaN values.
691 return
693 if sources_is_astropy:
694 psf_cat = sources[psf_mask]
695 else:
696 psf_cat = sources[psf_mask].copy(deep=True)
698 starXX = psf_cat[self.config.starShape + '_xx']
699 starYY = psf_cat[self.config.starShape + '_yy']
700 starXY = psf_cat[self.config.starShape + '_xy']
701 psfXX = psf_cat[self.config.psfShape + '_xx']
702 psfYY = psf_cat[self.config.psfShape + '_yy']
703 psfXY = psf_cat[self.config.psfShape + '_xy']
705 starE1 = (starXX - starYY)/(starXX + starYY)
706 starE2 = 2*starXY/(starXX + starYY)
708 psfE1 = (psfXX - psfYY)/(psfXX + psfYY)
709 psfE2 = 2*psfXY/(psfXX + psfYY)
711 e1Residuals = starE1 - psfE1
712 e2Residuals = starE2 - psfE2
714 # Comp TEx
715 ra = psf_cat["coord_ra"].to(units.deg)
716 dec = psf_cat["coord_dec"].to(units.deg)
718 psfTEx = {
719 "TE1": self.psfTE1,
720 "TE2": self.psfTE2,
721 }
723 gatherE12Stat = {
724 "TE1": {"E1": np.nan, "E2": np.nan, "Ex": np.nan, },
725 "TE2": {"E1": np.nan, "E2": np.nan, "Ex": np.nan, },
726 }
728 isNotNan = np.array([True] * len(ra))
729 isNotNan &= np.isfinite(ra)
730 isNotNan &= np.isfinite(dec)
731 isNotNan &= np.isfinite(e1Residuals)
732 isNotNan &= np.isfinite(e2Residuals)
734 if np.sum(isNotNan) >= 2:
735 # TE1 and TE2 computation, over visit.
736 for TEX in ["TE1", "TE2"]:
738 output = psfTEx[TEX].run(
739 e1Residuals[isNotNan], e2Residuals[isNotNan],
740 ra[isNotNan], dec[isNotNan],
741 units="degree",
742 )
744 gatherE12Stat[TEX]["E1"] = output.metric_E1
745 gatherE12Stat[TEX]["E2"] = output.metric_E2
746 gatherE12Stat[TEX]["Ex"] = output.metric_Ex
748 # TE3 and TE4 loop over detector and then median on visit.
750 psfTEx = {
751 "TE3": self.psfTE3,
752 "TE4": self.psfTE4,
753 }
755 gatherE34Stat = {
756 "TE3": {"E1": [], "E2": [], "Ex": [], },
757 "TE4": {"E1": [], "E2": [], "Ex": [], },
758 }
759 # calibrateImage run at detector level,
760 # need to wait second run of PSF to run this.
761 if "detector" in psf_cat.colnames:
762 detectorIds = list(set(psf_cat["detector"]))
763 for TEX in ["TE3", "TE4"]:
764 for ccdId in detectorIds:
765 isccdId = (ccdId == psf_cat["detector"])
766 mask = (isccdId & isNotNan)
767 if np.sum(mask) >= 2:
768 output = psfTEx[TEX].run(
769 e1Residuals[mask], e2Residuals[mask],
770 ra[mask], dec[mask],
771 units="degree",
772 )
773 gatherE34Stat[TEX]["E1"].append(output.metric_E1)
774 gatherE34Stat[TEX]["E2"].append(output.metric_E2)
775 gatherE34Stat[TEX]["Ex"].append(output.metric_Ex)
777 self.psfTE1e1 = summary.psfTE1e1 = gatherE12Stat["TE1"]["E1"]
778 self.psfTE1e2 = summary.psfTE1e2 = gatherE12Stat["TE1"]["E2"]
779 self.psfTE1ex = summary.psfTE1ex = gatherE12Stat["TE1"]["Ex"]
780 self.psfTE2e1 = summary.psfTE2e1 = gatherE12Stat["TE2"]["E1"]
781 self.psfTE2e2 = summary.psfTE2e2 = gatherE12Stat["TE2"]["E2"]
782 self.psfTE2ex = summary.psfTE2ex = gatherE12Stat["TE2"]["Ex"]
783 self.psfTE3e1 = summary.psfTE3e1 = np.nanmedian(gatherE34Stat["TE3"]["E1"])
784 self.psfTE3e2 = summary.psfTE3e2 = np.nanmedian(gatherE34Stat["TE3"]["E2"])
785 self.psfTE3ex = summary.psfTE3ex = np.nanmedian(gatherE34Stat["TE3"]["Ex"])
786 self.psfTE4e1 = summary.psfTE4e1 = np.nanmedian(gatherE34Stat["TE4"]["E1"])
787 self.psfTE4e2 = summary.psfTE4e2 = np.nanmedian(gatherE34Stat["TE4"]["E2"])
788 self.psfTE4ex = summary.psfTE4ex = np.nanmedian(gatherE34Stat["TE4"]["Ex"])
790 def update_wcs_stats(self, summary, wcs, bbox, visitInfo):
791 """Compute all summary-statistic fields that depend on the WCS model.
793 Parameters
794 ----------
795 summary : `lsst.afw.image.ExposureSummaryStats`
796 Summary object to update in-place.
797 wcs : `lsst.afw.geom.SkyWcs` or `None`
798 Astrometric calibration model. If `None`, all fields that depend
799 on the WCS will be reset (generally to NaN).
800 bbox : `lsst.geom.Box2I`
801 Bounding box of the image for which summary stats are being
802 computed.
803 visitInfo : `lsst.afw.image.VisitInfo`
804 Observation information used in together with ``wcs`` to compute
805 the zenith distance.
806 """
807 nan = float("nan")
808 summary.raCorners = [nan]*4
809 summary.decCorners = [nan]*4
810 summary.ra = nan
811 summary.dec = nan
812 summary.pixelScale = nan
813 summary.zenithDistance = nan
815 if wcs is None:
816 return
818 sph_pts = wcs.pixelToSky(geom.Box2D(bbox).getCorners())
819 summary.raCorners = [float(sph.getRa().asDegrees()) for sph in sph_pts]
820 summary.decCorners = [float(sph.getDec().asDegrees()) for sph in sph_pts]
822 sph_pt = wcs.pixelToSky(bbox.getCenter())
823 summary.ra = sph_pt.getRa().asDegrees()
824 summary.dec = sph_pt.getDec().asDegrees()
825 summary.pixelScale = wcs.getPixelScale(bbox.getCenter()).asArcseconds()
827 date = visitInfo.getDate()
829 if date.isValid(): 829 ↛ exitline 829 didn't return from function 'update_wcs_stats' because the condition on line 829 was always true
830 # We compute the zenithDistance at the center of the detector
831 # rather than use the boresight value available via the visitInfo,
832 # because the zenithDistance may vary significantly over a large
833 # field of view.
834 observatory = visitInfo.getObservatory()
835 loc = EarthLocation(lat=observatory.getLatitude().asDegrees()*units.deg,
836 lon=observatory.getLongitude().asDegrees()*units.deg,
837 height=observatory.getElevation()*units.m)
838 obstime = Time(visitInfo.getDate().get(system=DateTime.MJD),
839 location=loc, format="mjd")
840 coord = SkyCoord(
841 summary.ra*units.degree,
842 summary.dec*units.degree,
843 obstime=obstime,
844 location=loc,
845 )
846 with warnings.catch_warnings():
847 warnings.simplefilter("ignore")
848 altaz = coord.transform_to(AltAz)
850 summary.zenithDistance = float(90.0 - altaz.alt.degree)
852 def update_photo_calib_stats(self, summary, photo_calib):
853 """Compute all summary-statistic fields that depend on the photometric
854 calibration model.
856 Parameters
857 ----------
858 summary : `lsst.afw.image.ExposureSummaryStats`
859 Summary object to update in-place.
860 photo_calib : `lsst.afw.image.PhotoCalib` or `None`
861 Photometric calibration model. If `None`, all fields that depend
862 on the photometric calibration will be reset (generally to NaN).
863 """
864 if photo_calib is not None:
865 summary.zeroPoint = float(2.5*np.log10(photo_calib.getInstFluxAtZeroMagnitude()))
866 else:
867 summary.zeroPoint = float("nan")
869 def update_background_stats(self, summary, background):
870 """Compute summary-statistic fields that depend only on the
871 background model.
873 Parameters
874 ----------
875 summary : `lsst.afw.image.ExposureSummaryStats`
876 Summary object to update in-place.
877 background : `lsst.afw.math.BackgroundList` or `None`
878 Background model. If `None`, all fields that depend on the
879 background will be reset (generally to NaN).
881 Notes
882 -----
883 This does not include fields that depend on the background-subtracted
884 masked image; when the background changes, it should generally be
885 applied to the image and `update_masked_image_stats` should be called
886 as well.
887 """
888 if background is not None: 888 ↛ 893line 888 didn't jump to line 893 because the condition on line 888 was always true
889 bgStats = (bg[0].getStatsImage().getImage().array
890 for bg in background)
891 summary.skyBg = float(sum(np.median(bg[np.isfinite(bg)]) for bg in bgStats))
892 else:
893 summary.skyBg = float("nan")
895 def update_masked_image_stats(self, summary, masked_image):
896 """Compute summary-statistic fields that depend on the masked image
897 itself.
899 Parameters
900 ----------
901 summary : `lsst.afw.image.ExposureSummaryStats`
902 Summary object to update in-place.
903 masked_image : `lsst.afw.image.MaskedImage` or `None`
904 Masked image. If `None`, all fields that depend
905 on the masked image will be reset (generally to NaN).
906 """
907 nan = float("nan")
908 if masked_image is None: 908 ↛ 909line 908 didn't jump to line 909 because the condition on line 908 was never true
909 summary.skyNoise = nan
910 summary.meanVar = nan
911 return
912 statsCtrl = afwMath.StatisticsControl()
913 statsCtrl.setNumSigmaClip(self.config.sigmaClip)
914 statsCtrl.setNumIter(self.config.clipIter)
915 statsCtrl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.badMaskPlanes))
916 statsCtrl.setNanSafe(True)
918 statObj = afwMath.makeStatistics(masked_image, afwMath.STDEVCLIP, statsCtrl)
919 skyNoise, _ = statObj.getResult(afwMath.STDEVCLIP)
920 summary.skyNoise = skyNoise
922 statObj = afwMath.makeStatistics(masked_image.variance, masked_image.mask, afwMath.MEANCLIP,
923 statsCtrl)
924 meanVar, _ = statObj.getResult(afwMath.MEANCLIP)
925 summary.meanVar = meanVar
927 def update_effective_time_stats(self, summary, exposure):
928 """Compute effective exposure time statistics to estimate depth.
930 The effective exposure time is the equivalent shutter open
931 time that would be needed under nominal conditions to give the
932 same signal-to-noise for a point source as what is achieved by
933 the observation of interest. This metric combines measurements
934 of the point-spread function, the sky brightness, and the
935 transparency. It assumes that the observation is
936 sky-background dominated.
938 .. _teff_definitions:
940 The effective exposure time and its subcomponents are defined in [1]_.
942 References
943 ----------
945 .. [1] Neilsen, E.H., Bernstein, G., Gruendl, R., and Kent, S. (2016).
946 Limiting Magnitude, \tau, teff, and Image Quality in DES Year 1
947 https://www.osti.gov/biblio/1250877/
949 Parameters
950 ----------
951 summary : `lsst.afw.image.ExposureSummaryStats`
952 Summary object to update in-place.
953 exposure : `lsst.afw.image.ExposureF`
954 Exposure to grab band and exposure time metadata
956 """
957 nan = float("nan")
958 summary.effTime = nan
959 summary.effTimePsfSigmaScale = nan
960 summary.effTimeSkyBgScale = nan
961 summary.effTimeZeroPointScale = nan
963 exposureTime = exposure.getInfo().getVisitInfo().getExposureTime()
964 filterLabel = exposure.getFilter()
965 if (filterLabel is None) or (not filterLabel.hasBandLabel): 965 ↛ 966line 965 didn't jump to line 966 because the condition on line 965 was never true
966 band = None
967 else:
968 band = filterLabel.bandLabel
970 if band is None: 970 ↛ 971line 970 didn't jump to line 971 because the condition on line 970 was never true
971 self.log.warning("No band associated with exposure; effTime not calculated.")
972 return
974 # PSF component
975 if np.isnan(summary.psfSigma): 975 ↛ 976line 975 didn't jump to line 976 because the condition on line 975 was never true
976 self.log.debug("PSF sigma is NaN")
977 f_eff = nan
978 elif band not in self.config.fiducialPsfSigma:
979 self.log.debug(f"Fiducial PSF value not found for {band}")
980 f_eff = nan
981 else:
982 fiducialPsfSigma = self.config.fiducialPsfSigma[band]
983 f_eff = (summary.psfSigma / fiducialPsfSigma)**-2
985 # Transparency component
986 # Note: Assumes that the zeropoint includes the exposure time
987 if np.isnan(summary.zeroPoint):
988 self.log.debug("Zero point is NaN")
989 c_eff = nan
990 elif band not in self.config.fiducialZeroPoint:
991 self.log.debug(f"Fiducial zero point value not found for {band}")
992 c_eff = nan
993 else:
994 fiducialZeroPoint = self.config.fiducialZeroPoint[band]
995 zeroPointDiff = fiducialZeroPoint - (summary.zeroPoint - 2.5*np.log10(exposureTime))
996 c_eff = min(10**(-2.0*(zeroPointDiff)/2.5), self.config.maxEffectiveTransparency)
998 # Sky brightness component (convert to cts/s)
999 if np.isnan(summary.skyBg): 999 ↛ 1000line 999 didn't jump to line 1000 because the condition on line 999 was never true
1000 self.log.debug("Sky background is NaN")
1001 b_eff = nan
1002 elif band not in self.config.fiducialSkyBackground:
1003 self.log.debug(f"Fiducial sky background value not found for {band}")
1004 b_eff = nan
1005 else:
1006 fiducialSkyBackground = self.config.fiducialSkyBackground[band]
1007 b_eff = fiducialSkyBackground/(summary.skyBg/exposureTime)
1009 # Old effective exposure time scale factor
1010 # t_eff = f_eff * c_eff * b_eff
1012 # New effective exposure time (seconds) from magnitude limit
1013 if band not in self.config.fiducialMagLim:
1014 self.log.debug(f"Fiducial magnitude limit not found for {band}")
1015 effectiveTime = nan
1016 elif band not in self.config.fiducialExpTime: 1016 ↛ 1017line 1016 didn't jump to line 1017 because the condition on line 1016 was never true
1017 self.log.debug(f"Fiducial exposure time not found for {band}")
1018 effectiveTime = nan
1019 else:
1020 effectiveTime = compute_effective_time(summary.magLim,
1021 self.config.fiducialMagLim[band],
1022 self.config.fiducialExpTime[band])
1024 # Output quantities
1025 summary.effTime = float(effectiveTime)
1026 summary.effTimePsfSigmaScale = float(f_eff)
1027 summary.effTimeSkyBgScale = float(b_eff)
1028 summary.effTimeZeroPointScale = float(c_eff)
1030 def update_magnitude_limit_stats(self, summary, exposure):
1031 """Compute a summary statistic for the point-source magnitude limit at
1032 a given signal-to-noise ratio using exposure-level metadata.
1034 The magnitude limit is calculated at a given SNR from the PSF,
1035 sky, zeropoint, and readnoise in accordance with SMTN-002 [1]_,
1036 LSE-40 [2]_, and DMTN-296 [3]_.
1038 References
1039 ----------
1041 .. [1] "Calculating LSST limiting magnitudes and SNR" (SMTN-002)
1042 .. [2] "Photon Rates and SNR Calculations" (LSE-40)
1043 .. [3] "Calculations of Image and Catalog Depth" (DMTN-296)
1046 Parameters
1047 ----------
1048 summary : `lsst.afw.image.ExposureSummaryStats`
1049 Summary object to update in-place.
1050 exposure : `lsst.afw.image.ExposureF`
1051 Exposure to grab band and exposure time metadata
1052 """
1053 if exposure.getDetector() is None: 1053 ↛ 1054line 1053 didn't jump to line 1054 because the condition on line 1053 was never true
1054 summary.magLim = float("nan")
1055 return
1057 # Calculate the average readnoise [e-]
1058 readNoiseList = list(ipIsr.getExposureReadNoises(exposure).values())
1059 readNoise = np.nanmean(readNoiseList)
1060 if np.isnan(readNoise): 1060 ↛ 1061line 1060 didn't jump to line 1061 because the condition on line 1060 was never true
1061 readNoise = 0.0
1062 self.log.warning("Read noise set to NaN! Setting readNoise to 0.0 to proceed.")
1064 # Calculate the average gain [e-/ADU]
1065 gainList = list(ipIsr.getExposureGains(exposure).values())
1066 gain = np.nanmean(gainList)
1067 if np.isnan(gain): 1067 ↛ 1068line 1067 didn't jump to line 1068 because the condition on line 1067 was never true
1068 self.log.warning("Gain set to NaN! Setting magLim to NaN.")
1069 summary.magLim = float("nan")
1070 return
1072 # Get the image units (default to "adu" if metadata key absent)
1073 md = exposure.getMetadata()
1074 if md.get("LSST ISR UNITS", "adu") == "electron":
1075 gain = 1.0
1077 # Convert readNoise to image units
1078 readNoise /= gain
1080 # Calculate the limiting magnitude.
1081 # Note 1: Assumes that the image and readnoise have the same units
1082 # Note 2: Assumes that the zeropoint includes the exposure time
1083 magLim = compute_magnitude_limit(summary.psfArea, summary.skyBg,
1084 summary.zeroPoint, readNoise,
1085 gain, self.config.magLimSnr)
1087 summary.magLim = float(magLim)
1090def maximum_nearest_psf_distance(
1091 image_mask,
1092 psf_cat,
1093 sampling=8,
1094 bad_mask_bits=["BAD", "CR", "INTRP", "SAT", "SUSPECT", "NO_DATA", "EDGE"],
1095):
1096 """Compute the maximum distance of an unmasked pixel to its nearest PSF.
1098 Parameters
1099 ----------
1100 image_mask : `lsst.afw.image.Mask`
1101 The mask plane associated with the exposure.
1102 psf_cat : `lsst.afw.table.SourceCatalog` or `astropy.table.Table`
1103 Catalog containing only the stars used in the PSF modeling.
1104 sampling : `int`
1105 Sampling rate in each dimension to create the grid of points on which
1106 to evaluate the distance to the nearest PSF star. The tradeoff is
1107 between adequate sampling versus speed.
1108 bad_mask_bits : `list` [`str`]
1109 Mask bits required to be absent for a pixel to be considered
1110 "unmasked".
1112 Returns
1113 -------
1114 max_dist_to_nearest_psf : `float`
1115 The maximum distance (in pixels) of an unmasked pixel to its nearest
1116 PSF model star.
1117 """
1118 mask_arr = image_mask.array[::sampling, ::sampling]
1119 bitmask = image_mask.getPlaneBitMask(bad_mask_bits)
1120 good = ((mask_arr & bitmask) == 0)
1122 x = np.arange(good.shape[1]) * sampling
1123 y = np.arange(good.shape[0]) * sampling
1124 xx, yy = np.meshgrid(x, y)
1126 dist_to_nearest_psf = np.full(good.shape, np.inf)
1127 for psf in psf_cat:
1128 x_psf = psf["slot_Centroid_x"]
1129 y_psf = psf["slot_Centroid_y"]
1130 dist_to_nearest_psf = np.minimum(dist_to_nearest_psf, np.hypot(xx - x_psf, yy - y_psf))
1131 unmasked_dists = dist_to_nearest_psf * good
1132 max_dist_to_nearest_psf = np.max(unmasked_dists)
1134 return max_dist_to_nearest_psf
1137def compute_psf_image_deltas(
1138 image_mask,
1139 image_psf,
1140 sampling=96,
1141 ap_radius_pix=3.0,
1142 bad_mask_bits=["BAD", "CR", "INTRP", "SAT", "SUSPECT", "NO_DATA", "EDGE"],
1143):
1144 """Compute the delta between the maximum and minimum model PSF trace radius
1145 values evaluated on a grid of points lying in the unmasked region of the
1146 image.
1148 Parameters
1149 ----------
1150 image_mask : `lsst.afw.image.Mask`
1151 The mask plane associated with the exposure.
1152 image_psf : `lsst.afw.detection.Psf`
1153 The PSF model associated with the exposure.
1154 sampling : `int`, optional
1155 Sampling rate in each dimension to create the grid of points at which
1156 to evaluate ``image_psf``s trace radius value. The tradeoff is between
1157 adequate sampling versus speed.
1158 ap_radius_pix : `float`, optional
1159 Radius in pixels of the aperture on which to measure the flux of the
1160 PSF model.
1161 bad_mask_bits : `list` [`str`], optional
1162 Mask bits required to be absent for a pixel to be considered
1163 "unmasked".
1165 Returns
1166 -------
1167 psf_trace_radius_delta, psf_ap_flux_delta : `float`
1168 The delta (in pixels) between the maximum and minimum model PSF trace
1169 radius values and the PSF aperture fluxes (with aperture radius of
1170 max(2, 3*psfSigma)) evaluated on the x,y-grid subsampled on the
1171 unmasked detector pixels by a factor of ``sampling``. If both the
1172 model PSF trace radius value and aperture flux value on the grid
1173 evaluate to NaN, then NaNs are returned immediately.
1174 """
1175 psf_trace_radius_list = []
1176 psf_ap_flux_list = []
1177 mask_arr = image_mask.array[::sampling, ::sampling]
1178 bitmask = image_mask.getPlaneBitMask(bad_mask_bits)
1179 good = ((mask_arr & bitmask) == 0)
1181 x = np.arange(good.shape[1]) * sampling
1182 y = np.arange(good.shape[0]) * sampling
1183 xx, yy = np.meshgrid(x, y)
1185 for x_mesh, y_mesh, good_mesh in zip(xx, yy, good):
1186 for x_point, y_point, is_good in zip(x_mesh, y_mesh, good_mesh):
1187 if is_good:
1188 position = geom.Point2D(x_point, y_point)
1189 psf_trace_radius = image_psf.computeShape(position).getTraceRadius()
1190 psf_ap_flux = image_psf.computeApertureFlux(ap_radius_pix, position)
1191 if ~np.isfinite(psf_trace_radius) and ~np.isfinite(psf_ap_flux): 1191 ↛ 1192line 1191 didn't jump to line 1192 because the condition on line 1191 was never true
1192 return float("nan"), float("nan")
1193 psf_trace_radius_list.append(psf_trace_radius)
1194 psf_ap_flux_list.append(psf_ap_flux)
1196 psf_trace_radius_delta = np.max(psf_trace_radius_list) - np.min(psf_trace_radius_list)
1197 if np.any(np.asarray(psf_ap_flux_list) < 0.0): # Consider any -ve flux entry as "bad". 1197 ↛ 1198line 1197 didn't jump to line 1198 because the condition on line 1197 was never true
1198 psf_ap_flux_delta = float("nan")
1199 else:
1200 psf_ap_flux_delta = np.max(psf_ap_flux_list) - np.min(psf_ap_flux_list)
1202 return psf_trace_radius_delta, psf_ap_flux_delta
1205def compute_ap_corr_sigma_scaled_delta(
1206 image_mask,
1207 image_ap_corr_field,
1208 psfSigma,
1209 sampling=96,
1210 bad_mask_bits=["BAD", "CR", "INTRP", "SAT", "SUSPECT", "NO_DATA", "EDGE"],
1211):
1212 """Compute the delta between the maximum and minimum aperture correction
1213 values scaled (divided) by ``psfSigma`` for the given field representation,
1214 ``image_ap_corr_field`` evaluated on a grid of points lying in the
1215 unmasked region of the image.
1217 Parameters
1218 ----------
1219 image_mask : `lsst.afw.image.Mask`
1220 The mask plane associated with the exposure.
1221 image_ap_corr_field : `lsst.afw.math.ChebyshevBoundedField`
1222 The ChebyshevBoundedField representation of the aperture correction
1223 of interest for the exposure.
1224 psfSigma : `float`
1225 The PSF model second-moments determinant radius (center of chip)
1226 in pixels.
1227 sampling : `int`, optional
1228 Sampling rate in each dimension to create the grid of points at which
1229 to evaluate ``image_psf``s trace radius value. The tradeoff is between
1230 adequate sampling versus speed.
1231 bad_mask_bits : `list` [`str`], optional
1232 Mask bits required to be absent for a pixel to be considered
1233 "unmasked".
1235 Returns
1236 -------
1237 ap_corr_sigma_scaled_delta : `float`
1238 The delta between the maximum and minimum of the (multiplicative)
1239 aperture correction values scaled (divided) by ``psfSigma`` evaluated
1240 on the x,y-grid subsampled on the unmasked detector pixels by a factor
1241 of ``sampling``. If the aperture correction evaluates to NaN on any
1242 of the grid points, this is set to NaN.
1243 """
1244 mask_arr = image_mask.array[::sampling, ::sampling]
1245 bitmask = image_mask.getPlaneBitMask(bad_mask_bits)
1246 good = ((mask_arr & bitmask) == 0)
1248 x = np.arange(good.shape[1], dtype=np.float64) * sampling
1249 y = np.arange(good.shape[0], dtype=np.float64) * sampling
1250 xx, yy = np.meshgrid(x, y)
1252 ap_corr = image_ap_corr_field.evaluate(xx.ravel(), yy.ravel()).reshape(xx.shape)
1253 ap_corr_good = ap_corr[good]
1254 if ~np.isfinite(ap_corr_good).all(): 1254 ↛ 1255line 1254 didn't jump to line 1255 because the condition on line 1254 was never true
1255 ap_corr_sigma_scaled_delta = float("nan")
1256 else:
1257 ap_corr_sigma_scaled_delta = (np.max(ap_corr_good) - np.min(ap_corr_good))/psfSigma
1259 return ap_corr_sigma_scaled_delta
1262def compute_magnitude_limit(
1263 psfArea,
1264 skyBg,
1265 zeroPoint,
1266 readNoise,
1267 gain,
1268 snr
1269):
1270 """Compute the expected point-source magnitude limit at a given
1271 signal-to-noise ratio given the exposure-level metadata. Based on
1272 the signal-to-noise formula provided in SMTN-002 (see LSE-40 for
1273 more details on the calculation).
1275 SNR = C / sqrt( C/g + (B/g + sigma_inst**2) * neff )
1277 where C is the counts from the source, B is counts from the (sky)
1278 background, sigma_inst is the instrumental (read) noise, neff is
1279 the effective size of the PSF, and g is the gain in e-/ADU. Note
1280 that input values of ``skyBg``, ``zeroPoint``, and ``readNoise``
1281 should all consistently be in electrons or ADU.
1283 Parameters
1284 ----------
1285 psfArea : `float`
1286 The effective area of the PSF [pix].
1287 skyBg : `float`
1288 The sky background counts for the exposure [ADU or e-].
1289 zeroPoint : `float`
1290 The zeropoint (includes exposure time) [ADU or e-].
1291 readNoise : `float`
1292 The instrumental read noise for the exposure [ADU or e-].
1293 gain : `float`
1294 The instrumental gain for the exposure [e-/ADU]. The gain should
1295 be 1.0 if the skyBg, zeroPoint, and readNoise are in e-.
1296 snr : `float`
1297 Signal-to-noise ratio at which magnitude limit is calculated.
1299 Returns
1300 -------
1301 magnitude_limit : `float`
1302 The expected magnitude limit at the given signal to noise.
1303 """
1304 # Effective number of pixels within the PSF calculated directly
1305 # from the PSF model.
1306 neff = psfArea
1308 # Instrumental noise (read noise only)
1309 sigma_inst = readNoise
1311 # Total background counts derived from Eq. 45 of LSE-40
1312 background = (skyBg/gain + sigma_inst**2) * neff
1313 # Source counts to achieve the desired SNR (from quadratic formula)
1314 # Note typo in gain exponent in Eq. 45 of LSE-40 (DM-53234)
1315 source = (snr**2)/(2*gain) + np.sqrt((snr**4)/(4*gain**2) + snr**2 * background)
1317 # Convert to a magnitude using the zeropoint.
1318 # Note: Zeropoint currently includes exposure time
1319 magnitude_limit = -2.5*np.log10(source) + zeroPoint
1321 return magnitude_limit
1324def psf_sigma_to_psf_area(psfSigma, pixelScale):
1325 """Convert psfSigma [pixels] to an approximate psfArea [pixel^2] as defined in LSE-40.
1326 This is the same implementation followed by SMTN-002 to estimate SNR=5 magnitude limits.
1328 Parameters
1329 ----------
1330 psfSigma : `float`
1331 The PSF sigma value [pix].
1332 pixelScale : `float`
1333 The pixel scale [arcsec/pix].
1335 Returns
1336 -------
1337 psf_area : `float`
1338 Approximation of the PSF area [pix^2]
1339 """
1340 # Follow SMTN-002 to convert to geometric and effective FWHM
1341 fwhm_geom = psfSigma * 2.355 * pixelScale
1342 fwhm_eff = (fwhm_geom - 0.052) / 0.822
1343 # Area prefactor comes from LSE-40
1344 psf_area = 2.266 * (fwhm_eff / pixelScale)**2
1345 return psf_area
1348def compute_effective_time(
1349 magLim,
1350 fiducialMagLim,
1351 fiducialExpTime
1352):
1353 """Compute the effective exposure time from m5 following the prescription described in SMTN-296.
1355 teff = 10**(0.8 * (magLim - fiducialMagLim) ) * fiducialExpTime
1357 where `magLim` is the magnitude limit, `fiducialMagLim` is the fiducial magnitude limit calculated from
1358 the fiducial values of the ``psfArea``, ``skyBg``, ``zeroPoint``, and ``readNoise``, and
1359 `fiducialExpTime` is the fiducial exposure time (s).
1361 Parameters
1362 ----------
1363 magLim : `float`
1364 The measured magnitude limit [mag].
1365 fiducialMagLim : `float`
1366 The fiducial magnitude limit [mag].
1367 fiducialExpTime : `float`
1368 The fiducial exposure time [s].
1370 Returns
1371 -------
1372 effectiveTime : `float`
1373 The effective exposure time.
1374 """
1375 effectiveTime = 10**(0.8 * (magLim - fiducialMagLim)) * fiducialExpTime
1376 return effectiveTime