Coverage for python / lsst / pipe / tasks / peekExposure.py: 23%
418 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-05 18:37 +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/>.
21from __future__ import annotations
23__all__ = [
24 "PeekExposureTaskConfig",
25 "PeekExposureTask",
26]
28from typing import Any
30import astropy
31import numpy as np
32import numpy.typing as npt
34import lsst.afw.display as afwDisplay
35import lsst.afw.geom as afwGeom
36import lsst.afw.image as afwImage
37import lsst.afw.math as afwMath
38import lsst.afw.table as afwTable
39import lsst.daf.base as dafBase
40import lsst.geom as geom
41import lsst.pex.config as pexConfig
42import lsst.pipe.base as pipeBase
43from lsst.afw.detection import Psf
44from lsst.afw.geom.ellipses import Quadrupole
45from lsst.afw.image import ImageD
46from lsst.afw.table import SourceTable
47from lsst.geom import Box2I, Extent2I, LinearTransform, Point2D, Point2I, SpherePoint, arcseconds, degrees
48from lsst.meas.algorithms import SourceDetectionTask, SubtractBackgroundTask
49from lsst.meas.algorithms.installGaussianPsf import InstallGaussianPsfTask
50from lsst.meas.base import IdGenerator, SingleFrameMeasurementTask
51from lsst.meas.extensions import shapeHSM # noqa: F401
53IDX_SENTINEL = -99999
55# Adding obs_lsst to the table file isn't possible, so hard code this because
56# importing it from obs_lsst isn't possible without a deferred import, which
57# isn't possible because of the tests and the build order.
58FILTER_DELIMITER = "~"
61def isDispersedExp(exp):
62 """Check if an exposure is dispersed.
64 Note this is copied from `atmospec.utils.isDispersedExp` to avoid a
65 circular import.
67 Parameters
68 ----------
69 exp : `lsst.afw.image.Exposure`
70 The exposure.
72 Returns
73 -------
74 isDispersed : `bool`
75 Whether it is a dispersed image or not.
76 """
77 filterFullName = exp.filter.physicalLabel
78 if FILTER_DELIMITER not in filterFullName:
79 raise RuntimeError(f"Error parsing filter name {filterFullName}")
80 filt, grating = filterFullName.split(FILTER_DELIMITER)
81 if grating.upper().startswith('EMPTY'):
82 return False
83 return True
86def _estimateMode(data: npt.NDArray[np.float64], frac: float = 0.5) -> float:
87 """Estimate the mode of a 1d distribution.
89 Finds the smallest interval containing the fraction ``frac`` of the data,
90 then takes the median of the values in that interval.
92 Parameters
93 ----------
94 data : array-like
95 1d array of data values
96 frac : float, optional
97 Fraction of data to include in the mode interval. Default is 0.5.
99 Returns
100 -------
101 mode : float
102 Estimated mode of the data.
103 """
105 data = data[np.isfinite(data)]
106 if len(data) == 0:
107 return np.nan
108 elif len(data) == 1:
109 return data[0]
111 data = np.sort(data)
112 interval = int(np.ceil(frac * len(data)))
113 spans = data[interval:] - data[:-interval]
114 start = np.argmin(spans)
115 return np.median(data[start: start + interval])
118def _bearingToUnitVector(
119 wcs: afwGeom.SkyWcs,
120 bearing: geom.Angle,
121 imagePoint: geom.Point2D,
122 skyPoint: geom.SpherePoint | None = None,
123) -> geom.Extent2D:
124 """Compute unit vector along given bearing at given point in the sky.
126 Parameters
127 ----------
128 wcs : `lsst.afw.geom.SkyWcs`
129 World Coordinate System of image.
130 bearing : `lsst.geom.Angle`
131 Bearing (angle North of East) at which to compute unit vector.
132 imagePoint : `lsst.geom.Point2D`
133 Point in the image.
134 skyPoint : `lsst.geom.SpherePoint`, optional
135 Point in the sky.
137 Returns
138 -------
139 unitVector : `lsst.geom.Extent2D`
140 Unit vector in the direction of bearing.
141 """
142 if skyPoint is None:
143 skyPoint = wcs.pixelToSky(imagePoint)
144 dpt = wcs.skyToPixel(skyPoint.offset(bearing, 1e-4 * degrees)) - imagePoint
145 return dpt / dpt.computeNorm()
148def roseVectors(wcs: afwGeom.SkyWcs, imagePoint: geom.Point2D, parAng: geom.Angle | None = None) -> dict:
149 """Compute unit vectors in the N/W and optionally alt/az directions.
151 Parameters
152 ----------
153 wcs : `lsst.afw.geom.SkyWcs`
154 World Coordinate System of image.
155 imagePoint : `lsst.geom.Point2D`
156 Point in the image
157 parAng : `lsst.geom.Angle`, optional
158 Parallactic angle (position angle of zenith measured East from North)
159 (default: None)
161 Returns
162 -------
163 unitVectors : `dict` of `lsst.geom.Extent2D`
164 Unit vectors in the N, W, alt, and az directions.
165 """
166 ncp = SpherePoint(0 * degrees, 90 * degrees) # North Celestial Pole
167 skyPoint = wcs.pixelToSky(imagePoint)
168 bearing = skyPoint.bearingTo(ncp)
170 out = dict()
171 out["N"] = _bearingToUnitVector(wcs, bearing, imagePoint, skyPoint=skyPoint)
172 out["W"] = _bearingToUnitVector(wcs, bearing + 90 * degrees, imagePoint, skyPoint=skyPoint)
174 if parAng is not None:
175 out["alt"] = _bearingToUnitVector(wcs, bearing - parAng, imagePoint, skyPoint=skyPoint)
176 out["az"] = _bearingToUnitVector(wcs, bearing - parAng + 90 * degrees, imagePoint, skyPoint=skyPoint)
178 return out
181def plotRose(
182 display: afwDisplay.Display,
183 wcs: afwGeom.SkyWcs,
184 imagePoint: geom.Point2D,
185 parAng: geom.Angle | None = None,
186 len: float = 50,
187) -> None:
188 """Display unit vectors along N/W and optionally alt/az directions.
190 Parameters
191 ----------
192 display : `lsst.afw.display.Display`
193 Display on which to render rose.
194 wcs : `lsst.afw.geom.SkyWcs`
195 World Coordinate System of image.
196 imagePoint : `lsst.geom.Point2D`
197 Point in the image at which to render rose.
198 parAng : `lsst.geom.Angle`, optional
199 Parallactic angle (position angle of zenith measured East from North)
200 (default: None)
201 len : `float`, optional
202 Length of the rose vectors (default: 50)
203 """
204 unitVectors = roseVectors(wcs, imagePoint, parAng=parAng)
205 colors = dict(N="r", W="r", alt="g", az="g")
206 for name, unitVector in unitVectors.items():
207 display.line([imagePoint, imagePoint + len * unitVector], ctype=colors[name])
208 display.dot(name, *(imagePoint + 1.6 * len * unitVector), ctype=colors[name])
211class DonutPsf(Psf):
212 def __init__(self, size: float, outerRad: float, innerRad: float):
213 Psf.__init__(self, isFixed=True)
214 self.size = size
215 self.outerRad = outerRad
216 self.innerRad = innerRad
217 self.dimensions = Extent2I(size, size)
219 def __deepcopy__(self, memo: Any | None = None) -> DonutPsf:
220 return DonutPsf(self.size, self.outerRad, self.innerRad)
222 def resized(self, width: float, height: float) -> DonutPsf:
223 assert width == height
224 return DonutPsf(width, self.outerRad, self.innerRad)
226 def _doComputeKernelImage(
227 self, position: Point2D | None = None, color: afwImage.Color | None = None
228 ) -> ImageD:
229 bbox = self.computeBBox(self.getAveragePosition())
230 img = ImageD(bbox, 0.0)
231 x, y = np.ogrid[bbox.minY: bbox.maxY + 1, bbox.minX: bbox.maxX + 1]
232 rsqr = x**2 + y**2
233 w = (rsqr < self.outerRad**2) & (rsqr > self.innerRad**2)
234 img.array[w] = 1.0
235 img.array /= np.sum(img.array)
236 return img
238 def _doComputeBBox(self, position: Point2D | None = None, color: afwImage.Color | None = None) -> Box2I:
239 return Box2I(Point2I(-self.dimensions / 2), self.dimensions)
241 def _doComputeShape(
242 self, position: Point2D | None = None, color: afwImage.Color | None = None
243 ) -> Quadrupole:
244 Ixx = self.outerRad**4 - self.innerRad**4
245 Ixx /= self.outerRad**2 - self.innerRad**2
246 return Quadrupole(Ixx, Ixx, 0.0)
248 def _doComputeApertureFlux(
249 self, radius: float, position: Point2D | None = None, color: afwImage.Color | None = None
250 ) -> float:
251 return 1 - np.exp(-0.5 * (radius / self.sigma) ** 2)
253 def __eq__(self, rhs: object) -> bool:
254 if isinstance(rhs, DonutPsf):
255 return self.size == rhs.size and self.outerRad == rhs.outerRad and self.innerRad == rhs.innerRad
256 return False
259class PeekTaskConfig(pexConfig.Config):
260 """Config class for the PeekTask."""
262 installPsf = pexConfig.ConfigurableField(
263 target=InstallGaussianPsfTask,
264 doc="Install a PSF model",
265 ) # type: ignore
266 doInstallPsf: pexConfig.Field[bool] = pexConfig.Field(
267 dtype=bool,
268 default=True,
269 doc="Install a PSF model?",
270 )
271 background = pexConfig.ConfigurableField(
272 target=SubtractBackgroundTask,
273 doc="Estimate background",
274 ) # type: ignore
275 detection = pexConfig.ConfigurableField(target=SourceDetectionTask, doc="Detect sources") # type: ignore
276 measurement = pexConfig.ConfigurableField(
277 target=SingleFrameMeasurementTask, doc="Measure sources"
278 ) # type: ignore
279 defaultBinSize: pexConfig.Field[int] = pexConfig.Field(
280 dtype=int,
281 default=1,
282 doc="Default binning factor for exposure (often overridden).",
283 )
285 def setDefaults(self) -> None:
286 super().setDefaults()
287 # Configure to be aggressively fast.
288 self.detection.thresholdValue = 5.0
289 self.detection.includeThresholdMultiplier = 10.0
290 self.detection.reEstimateBackground = False
291 self.detection.doTempLocalBackground = False
292 self.measurement.doReplaceWithNoise = False
293 self.detection.minPixels = 40
294 self.installPsf.fwhm = 5.0
295 self.installPsf.width = 21
296 # minimal set of measurements
297 self.measurement.plugins.names = [
298 "base_PixelFlags",
299 "base_SdssCentroid",
300 "ext_shapeHSM_HsmSourceMoments",
301 "base_GaussianFlux",
302 "base_PsfFlux",
303 "base_CircularApertureFlux",
304 ]
305 self.measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"
308class PeekTask(pipeBase.Task):
309 """Peek at exposure to quickly detect and measure both the brightest source
310 in the image, and a set of sources representative of the exposure's overall
311 image quality.
313 Optionally bins image and then:
314 - installs a simple PSF model
315 - measures and subtracts the background
316 - detects sources
317 - measures sources
319 Designed to be quick at the expense of primarily completeness, but also to
320 a lesser extent accuracy.
321 """
323 ConfigClass = PeekTaskConfig
324 config: PeekTaskConfig
325 installPsf: InstallGaussianPsfTask
326 background: SubtractBackgroundTask
327 detection: SourceDetectionTask
328 measurement: SingleFrameMeasurementTask
329 _DefaultName = "peek"
331 def __init__(self, schema: Any | None = None, **kwargs: Any):
332 super().__init__(**kwargs)
334 if schema is None:
335 schema = SourceTable.makeMinimalSchema()
336 self.schema = schema
338 self.makeSubtask("installPsf")
339 self.makeSubtask("background")
340 self.makeSubtask("detection", schema=self.schema)
341 self.algMetadata = dafBase.PropertyList()
342 self.makeSubtask("measurement", schema=self.schema, algMetadata=self.algMetadata)
344 def run(self, exposure: afwImage.Exposure, binSize: int | None = None) -> pipeBase.Struct:
345 """Peek at exposure.
347 Parameters
348 ----------
349 exposure : `lsst.afw.image.Exposure`
350 Exposure at which to peek.
351 binSize : `int`, optional
352 Binning factor for exposure. Default is None, which will use the
353 default binning factor from the config.
355 Returns
356 -------
357 result : `pipeBase.Struct`
358 Result of peeking.
359 Struct containing:
360 - sourceCat : `lsst.afw.table.SourceCatalog`
361 Source catalog from the binned exposure.
362 """
363 if binSize is None:
364 binSize = self.config.defaultBinSize
366 if binSize != 1:
367 mi = exposure.getMaskedImage()
368 binned = afwMath.binImage(mi, binSize)
369 exposure.setMaskedImage(binned)
371 if self.config.doInstallPsf:
372 self.installPsf.run(exposure=exposure)
374 self.background.run(exposure)
376 idGenerator = IdGenerator()
377 sourceIdFactory = idGenerator.make_table_id_factory()
378 table = SourceTable.make(self.schema, sourceIdFactory)
379 table.setMetadata(self.algMetadata)
380 sourceCat = self.detection.run(table=table, exposure=exposure, doSmooth=True).sources
382 self.measurement.run(measCat=sourceCat, exposure=exposure, exposureId=idGenerator.catalog_id)
384 return pipeBase.Struct(
385 sourceCat=sourceCat,
386 )
389class PeekDonutTaskConfig(pexConfig.Config):
390 """Config class for the PeekDonutTask."""
392 peek = pexConfig.ConfigurableField(
393 target=PeekTask,
394 doc="Peek configuration",
395 ) # type: ignore
396 resolution = pexConfig.Field(
397 dtype=float,
398 default=16.0,
399 doc="Target number of pixels for a binned donut",
400 ) # type: ignore
401 binSizeMax = pexConfig.Field(
402 dtype=int,
403 default=10,
404 doc="Maximum binning factor for donut mode",
405 ) # type: ignore
407 def setDefaults(self) -> None:
408 super().setDefaults()
409 # Donuts are big even when binned.
410 self.peek.installPsf.fwhm = 10.0
411 self.peek.installPsf.width = 31
412 # Use DonutPSF if not overridden
413 self.peek.doInstallPsf = False
416class PeekDonutTask(pipeBase.Task):
417 """Peek at a donut exposure.
419 The main modification for donuts is to aggressively bin the image to reduce
420 the size of sources (donuts) from ~100 pixels or more to ~10 pixels. This
421 greatly increases the speed and detection capabilities of PeekTask with
422 little loss of accuracy for centroids.
423 """
425 ConfigClass = PeekDonutTaskConfig
426 config: PeekDonutTaskConfig
427 peek: PeekTask
428 _DefaultName = "peekDonut"
430 def __init__(self, config: Any, **kwargs: Any):
431 super().__init__(config=config, **kwargs)
432 self.makeSubtask("peek")
434 def run(
435 self, exposure: afwImage.Exposure, donutDiameter: float, binSize: int | None = None
436 ) -> pipeBase.Struct:
437 """Peek at donut exposure.
439 Parameters
440 ----------
441 exposure : `lsst.afw.image.Exposure`
442 Exposure at which to peek.
443 donutDiameter : `float`
444 Donut diameter in pixels.
445 binSize : `int`, optional
446 Binning factor for exposure. Default is None, which will use the
447 resolution config value to determine the binSize.
449 Returns
450 -------
451 result : `pipeBase.Struct`
452 Result of donut peeking.
453 Struct containing:
454 - mode : `str`
455 Peek mode that was run.
456 - binSize : `int`
457 Binning factor used.
458 - binnedSourceCat : `lsst.afw.table.SourceCatalog`
459 Source catalog from the binned exposure.
460 """
461 if binSize is None:
462 binSize = int(
463 np.floor(
464 np.clip(
465 donutDiameter / self.config.resolution,
466 1,
467 self.config.binSizeMax,
468 )
469 )
470 )
471 binnedDonutDiameter = donutDiameter / binSize
472 psf = DonutPsf(
473 binnedDonutDiameter * 1.5, binnedDonutDiameter * 0.5, binnedDonutDiameter * 0.5 * 0.3525
474 )
476 # Note that SourceDetectionTask will convolve with a _Gaussian
477 # approximation to the PSF_ anyway, so we don't really need to be
478 # precise with the PSF unless this changes. PSFs that approach the
479 # size of the image, however, can cause problems with the detection
480 # convolution algorithm, so we limit the size.
481 sigma = psf.computeShape(psf.getAveragePosition()).getDeterminantRadius()
482 factor = 8 * sigma / (min(exposure.getDimensions()) / binSize)
484 if factor > 1:
485 psf = DonutPsf(
486 binnedDonutDiameter * 1.5 / factor,
487 binnedDonutDiameter * 0.5 / factor,
488 binnedDonutDiameter * 0.5 * 0.3525 / factor,
489 )
490 exposure.setPsf(psf)
492 peekResult = self.peek.run(exposure, binSize)
494 return pipeBase.Struct(
495 mode="donut",
496 binSize=binSize,
497 binnedSourceCat=peekResult.sourceCat,
498 )
500 def getGoodSources(self, binnedSourceCat: afwTable.SourceCatalog) -> np.ndarray:
501 """Perform any filtering on the source catalog.
503 Parameters
504 ----------
505 binnedSourceCat : `lsst.afw.table.SourceCatalog`
506 Source catalog from the binned exposure.
508 Returns
509 -------
510 goodSourceMask : `numpy.ndarray`
511 Boolean array indicating which sources are good.
512 """
513 # Perform any filtering on the source catalog
514 goodSourceMask = np.ones(len(binnedSourceCat), dtype=bool)
515 return goodSourceMask
518class PeekPhotoTaskConfig(pexConfig.Config):
519 """Config class for the PeekPhotoTask."""
521 peek = pexConfig.ConfigurableField(
522 target=PeekTask,
523 doc="Peek configuration",
524 ) # type: ignore
525 binSize: pexConfig.Field[int] = pexConfig.Field(
526 dtype=int,
527 default=2,
528 doc="Binning factor for exposure",
529 )
531 def setDefaults(self) -> None:
532 super().setDefaults()
533 # Use a lower detection threshold in photo mode to go a bit fainter.
534 self.peek.detection.includeThresholdMultiplier = 1.0
535 self.peek.detection.thresholdValue = 10.0
536 self.peek.detection.minPixels = 10
539class PeekPhotoTask(pipeBase.Task):
540 """Peek at a photo (imaging) exposure.
542 For photo mode, we keep a relatively small detection threshold value, so we
543 can detect faint sources to use for image quality assessment.
544 """
546 ConfigClass = PeekPhotoTaskConfig
547 config: PeekPhotoTaskConfig
548 peek: PeekTask
549 _DefaultName = "peekPhoto"
551 def __init__(self, config: Any, **kwargs: Any):
552 super().__init__(config=config, **kwargs)
553 self.makeSubtask("peek")
555 def run(self, exposure: afwImage.Exposure, binSize: int | None = None) -> pipeBase.Struct:
556 """Peek at donut exposure.
558 Parameters
559 ----------
560 exposure : `lsst.afw.image.Exposure`
561 Exposure at which to peek.
562 binSize : `int`, optional
563 Binning factor for exposure. Default is None, which will use the
564 binning factor from the config.
566 Returns
567 -------
568 result : `pipeBase.Struct`
569 Result of photo peeking.
570 Struct containing:
571 - mode : `str`
572 Peek mode that was run.
573 - binSize : `int`
574 Binning factor used.
575 - binnedSourceCat : `lsst.afw.table.SourceCatalog`
576 Source catalog from the binned exposure.
577 """
578 if binSize is None:
579 binSize = self.config.binSize
581 peekResult = self.peek.run(exposure, binSize)
583 return pipeBase.Struct(
584 mode="photo",
585 binSize=binSize,
586 binnedSourceCat=peekResult.sourceCat,
587 )
589 def getGoodSources(self, binnedSourceCat: afwTable.SourceCatalog) -> np.ndarray:
590 """Perform any filtering on the source catalog.
592 Parameters
593 ----------
594 binnedSourceCat : `lsst.afw.table.SourceCatalog`
595 Source catalog from the binned exposure.
597 Returns
598 -------
599 goodSourceMask : `numpy.ndarray`
600 Boolean array indicating which sources are good.
601 """
602 # Perform any filtering on the source catalog
603 goodSourceMask = np.ones(len(binnedSourceCat), dtype=bool)
604 return goodSourceMask
607class PeekSpecTaskConfig(pexConfig.Config):
608 """Config class for the PeekSpecTask."""
610 peek = pexConfig.ConfigurableField(
611 target=PeekTask,
612 doc="Peek configuration",
613 ) # type: ignore
614 binSize: pexConfig.Field[int] = pexConfig.Field(
615 dtype=int,
616 default=2,
617 doc="binning factor for exposure",
618 )
619 maxFootprintAspectRatio: pexConfig.Field[int] = pexConfig.Field(
620 dtype=float,
621 default=10.0,
622 doc="Maximum detection footprint aspect ratio to consider as 0th order" " (non-dispersed) light.",
623 )
625 def setDefaults(self) -> None:
626 super().setDefaults()
627 # Use bright threshold
628 self.peek.detection.includeThresholdMultiplier = 1.0
629 self.peek.detection.thresholdValue = 500.0
630 # Use a large radius aperture flux for spectra to better identify the
631 # brightest source, which for spectra often has a saturated core.
632 self.peek.measurement.slots.apFlux = "base_CircularApertureFlux_70_0"
633 # Also allow a larger distance to peak for centroiding in case there's
634 # a relatively large saturated region.
635 self.peek.measurement.plugins["base_SdssCentroid"].maxDistToPeak = 15.0
638class PeekSpecTask(pipeBase.Task):
639 """Peek at a spectroscopic exposure.
641 For spec mode, we dramatically increase the detection threshold to avoid
642 creating blends with the long spectra objects that appear in these images.
643 We also change the default aperture flux slot to a larger aperture, which
644 helps overcome challenges with lost flux in the interpolated cores of
645 saturated objects.
646 """
648 ConfigClass = PeekSpecTaskConfig
649 config: PeekSpecTaskConfig
650 peek: PeekTask
651 _DefaultName = "peekSpec"
653 def __init__(self, config: Any, **kwargs: Any):
654 super().__init__(config=config, **kwargs)
655 self.makeSubtask("peek")
657 def run(self, exposure: afwImage.Exposure, binSize: int | None = None) -> pipeBase.Struct:
658 """Peek at spectroscopic exposure.
660 Parameters
661 ----------
662 exposure : `lsst.afw.image.Exposure`
663 Exposure at which to peek.
664 binSize : `int`, optional
665 Binning factor for exposure. Default is None, which will use the
666 binning factor from the config.
668 Returns
669 -------
670 result : `pipeBase.Struct`
671 Result of spec peeking.
672 Struct containing:
673 - mode : `str`
674 Peek mode that was run.
675 - binSize : `int`
676 Binning factor used.
677 - binnedSourceCat : `lsst.afw.table.SourceCatalog`
678 Source catalog from the binned exposure.
679 """
680 if binSize is None:
681 binSize = self.config.binSize
683 peekResult = self.peek.run(exposure, binSize)
685 return pipeBase.Struct(
686 mode="spec",
687 binSize=binSize,
688 binnedSourceCat=peekResult.sourceCat,
689 )
691 def getGoodSources(self, binnedSourceCat: afwTable.SourceCatalog) -> np.ndarray:
692 """Perform any filtering on the source catalog.
694 Parameters
695 ----------
696 binnedSourceCat : `lsst.afw.table.SourceCatalog`
697 Source catalog from the binned exposure.
699 Returns
700 -------
701 goodSourceMask : `numpy.ndarray`
702 Boolean array indicating which sources are good.
703 """
704 # Perform any filtering on the source catalog
705 goodSourceMask = np.ones(len(binnedSourceCat), dtype=bool)
706 fpShapes = [src.getFootprint().getShape() for src in binnedSourceCat]
707 # Filter out likely spectrum detections
708 goodSourceMask &= np.array(
709 [sh.getIyy() < self.config.maxFootprintAspectRatio * sh.getIxx() for sh in fpShapes],
710 dtype=np.bool_,
711 )
712 return goodSourceMask
715class PeekExposureTaskConfig(pexConfig.Config):
716 """Config class for the PeekExposureTask."""
718 donutThreshold: pexConfig.Field[float] = pexConfig.Field(
719 dtype=float,
720 default=50.0,
721 doc="Size threshold in pixels for when to switch to donut mode.",
722 )
723 doPhotoFallback: pexConfig.Field[bool] = pexConfig.Field(
724 dtype=bool,
725 default=True,
726 doc="If True, fall back to photo mode if spec mode fails.",
727 )
728 donut = pexConfig.ConfigurableField(
729 target=PeekDonutTask,
730 doc="PeekDonut task",
731 ) # type: ignore
732 photo = pexConfig.ConfigurableField(
733 target=PeekPhotoTask,
734 doc="PeekPhoto task",
735 ) # type: ignore
736 spec = pexConfig.ConfigurableField(
737 target=PeekSpecTask,
738 doc="PeekSpec task",
739 ) # type: ignore
742class PeekExposureTask(pipeBase.Task):
743 """Peek at exposure to quickly detect and measure both the brightest
744 source in the image, and a set of sources representative of the
745 exposure's overall image quality.
747 Parameters
748 ----------
749 config : `lsst.summit.utils.peekExposure.PeekExposureTaskConfig`
750 Configuration for the task.
751 display : `lsst.afw.display.Display`, optional
752 For displaying the exposure and sources.
754 Notes
755 -----
756 The basic philosophy of PeekExposureTask is to:
757 1) Classify exposures based on metadata into 'donut', 'spec', or 'photo'.
758 2) Run PeekTask on the exposure through a wrapper with class specific
759 modifications.
760 3) Try only to branch in the code based on the metadata, and not on the
761 data itself. This avoids problematic discontinuities in measurements.
763 The main knobs we fiddle with based on the classification are:
764 - Detection threshold
765 - Minimum number of pixels for a detection
766 - Binning of the image
767 - Installed PSF size
768 """
770 ConfigClass = PeekExposureTaskConfig
771 config: PeekExposureTaskConfig
772 donut: PeekDonutTask
773 photo: PeekPhotoTask
774 spec: PeekSpecTask
775 _DefaultName = "peekExposureTask"
777 def __init__(self, config: Any, *, display: Any = None, **kwargs: Any):
778 super().__init__(config=config, **kwargs)
780 self.makeSubtask("donut")
781 self.makeSubtask("photo")
782 self.makeSubtask("spec")
784 self._display = display
786 def getDonutDiameter(self, exposure: afwImage.Exposure) -> float:
787 """Estimate donut diameter from exposure metadata.
789 Parameters
790 ----------
791 exposure : `lsst.afw.image.Exposure`
792 Exposure to estimate donut diameter for.
794 Returns
795 -------
796 donutDiameter : `float`
797 Estimated donut diameter in pixels.
798 """
799 visitInfo = exposure.getInfo().getVisitInfo()
800 focusZ = visitInfo.focusZ
801 instrumentLabel = visitInfo.instrumentLabel
803 match instrumentLabel:
804 case "LATISS":
805 focusZ *= 41 # magnification factor
806 fratio = 18.0
807 case "LSSTCam" | "LSSTComCam" | "LSSTComCamSim":
808 fratio = 1.234
809 case _:
810 raise ValueError(f"Unknown instrument label: {instrumentLabel}")
811 # AuxTel/ComCam/LSSTCam all have 10 micron pixels (= 10e-3 mm)
812 donutDiameter = abs(focusZ) / fratio / 10e-3
813 self.log.info(f"{focusZ=} mm")
814 self.log.info(f"donutDiameter = {donutDiameter} pixels")
815 return donutDiameter
817 def run(
818 self,
819 exposure: afwImage.Exposure,
820 *,
821 doDisplay: bool = False,
822 doDisplayIndices: bool = False,
823 mode: str = "auto",
824 binSize: int | None = None,
825 donutDiameter: float | None = None,
826 ) -> pipeBase.Struct:
827 """
828 Parameters
829 ----------
830 exposure : `lsst.afw.image.Exposure`
831 Exposure at which to peek.
832 doDisplay : `bool`, optional
833 Display the exposure and sources? Default False. (Requires
834 display to have been passed to task constructor)
835 doDisplayIndices : `bool`, optional
836 Display the source indices? Default False. (Requires display to
837 have been passed to task constructor)
838 mode : {'auto', 'donut', 'spec', 'photo'}, optional
839 Mode to run in. Default 'auto'.
840 binSize : `int`, optional
841 Binning factor for exposure. Default is None, which let's subtasks
842 control rebinning directly.
843 donutDiameter : `float`, optional
844 Donut diameter in pixels. Default is None, which will estimate the
845 donut diameter from the exposure metadata.
847 Returns
848 -------
849 result : `pipeBase.Struct`
850 Result of the peek.
851 Struct containing:
852 - mode : `str`
853 Peek mode that was run.
854 - binSize : `int`
855 Binning factor used.
856 - binnedSourceCat : `lsst.afw.table.SourceCatalog`
857 Source catalog from the binned exposure.
858 - table : `astropy.table.Table`
859 Curated source table in unbinned coordinates.
860 - brightestIdx : `int`
861 Index of brightest source in source catalog.
862 - brightestCentroid : `lsst.geom.Point2D`
863 Brightest source centroid in unbinned pixel coords.
864 - brightestPixelShape : `lsst.afw.geom.Quadrupole`
865 Brightest source shape in unbinned pixel coords.
866 - brightestEquatorialShape : `lsst.afw.geom.Quadrupole`
867 Brightest source shape in equitorial coordinates (arcsec).
868 - brightestAltAzShape : `lsst.afw.geom.Quadrupole`
869 Brightest source shape in alt/az coordinates (arcsec).
870 - psfPixelShape : `lsst.afw.geom.Quadrupole`
871 Estimated PSF shape in unbinned pixel coords.
872 - psfEquatorialShape : `lsst.afw.geom.Quadrupole`
873 Estimated PSF shape in equitorial coordinates (arcsec).
874 - psfAltAzShape : `lsst.afw.geom.Quadrupole`
875 Estimated PSF shape in alt/az coordinates (arcsec).
876 - pixelMedian : `float`
877 Median estimate of entire image.
878 - pixelMode : `float`
879 Mode estimate of entire image.
880 """
881 # Make a copy so the original image is unmodified.
882 exposure = exposure.clone()
883 try:
884 result = self._run(exposure, doDisplay, doDisplayIndices, mode, binSize, donutDiameter)
885 except Exception as e:
886 self.log.warning(f"Peek failed: {e}")
887 result = pipeBase.Struct(
888 mode="failed",
889 binSize=0,
890 binnedSourceCat=None,
891 table=None,
892 brightestIdx=0,
893 brightestCentroid=Point2D(np.nan, np.nan),
894 brightestPixelShape=Quadrupole(np.nan, np.nan, np.nan),
895 brightestEquatorialShape=Quadrupole(np.nan, np.nan, np.nan),
896 brightestAltAzShape=Quadrupole(np.nan, np.nan, np.nan),
897 psfPixelShape=Quadrupole(np.nan, np.nan, np.nan),
898 psfEquatorialShape=Quadrupole(np.nan, np.nan, np.nan),
899 psfAltAzShape=Quadrupole(np.nan, np.nan, np.nan),
900 pixelMedian=np.nan,
901 pixelMode=np.nan,
902 )
903 return result
905 def _run(
906 self,
907 exposure: afwImage.Exposure,
908 doDisplay: bool,
909 doDisplayIndices: bool,
910 mode: str,
911 binSize: int | None,
912 donutDiameter: float | None,
913 ) -> pipeBase.Struct:
914 """The actual run method, called by run()."""
915 # If image is ~large, then use a subsampling of the image for
916 # speedy median/mode estimates.
917 arr = exposure.getMaskedImage().getImage().array
918 sampling = 1
919 if arr.size > 250_000:
920 sampling = int(np.floor(np.sqrt(arr.size / 250_000)))
921 pixelMedian = np.nanmedian(arr[::sampling, ::sampling])
922 pixelMode = _estimateMode(arr[::sampling, ::sampling])
924 if donutDiameter is None:
925 donutDiameter = self.getDonutDiameter(exposure)
927 mode, binSize, binnedSourceCat = self.runPeek(exposure, mode, donutDiameter, binSize)
929 table = self.transformTable(binSize, binnedSourceCat)
931 match mode:
932 case "donut":
933 goodSourceMask = self.donut.getGoodSources(binnedSourceCat)
934 case "spec":
935 goodSourceMask = self.spec.getGoodSources(binnedSourceCat)
936 case "photo":
937 goodSourceMask = self.photo.getGoodSources(binnedSourceCat)
939 # prepare output variables
940 maxFluxIdx, brightCentroid, brightShape = self.getBrightest(binnedSourceCat, binSize, goodSourceMask)
941 psfShape = self.getPsfShape(binnedSourceCat, binSize, goodSourceMask)
943 equatorialShapes, altAzShapes = self.transformShapes([brightShape, psfShape], exposure, binSize)
945 if doDisplay:
946 self.updateDisplay(exposure, binSize, binnedSourceCat, maxFluxIdx, doDisplayIndices)
948 return pipeBase.Struct(
949 mode=mode,
950 binSize=binSize,
951 binnedSourceCat=binnedSourceCat,
952 table=table,
953 brightestIdx=maxFluxIdx,
954 brightestCentroid=brightCentroid,
955 brightestPixelShape=brightShape,
956 brightestEquatorialShape=equatorialShapes[0],
957 brightestAltAzShape=altAzShapes[0],
958 psfPixelShape=psfShape,
959 psfEquatorialShape=equatorialShapes[1],
960 psfAltAzShape=altAzShapes[1],
961 pixelMedian=pixelMedian,
962 pixelMode=pixelMode,
963 )
965 def runPeek(
966 self,
967 exposure: afwImage.Exposure,
968 mode: str,
969 donutDiameter: float,
970 binSize: int | None = None,
971 ) -> tuple[str, int, afwTable.SourceCatalog]:
972 """Classify exposure and run appropriate PeekTask wrapper.
974 Parameters
975 ----------
976 exposure : `lsst.afw.image.Exposure`
977 Exposure to peek.
978 mode : {'auto', 'donut', 'spec', 'photo'}
979 Mode to run in.
980 donutDiameter : `float`
981 Donut diameter in pixels.
982 binSize : `int`, optional
983 Binning factor for exposure. Default is None, which let's subtasks
984 control rebinning directly.
986 Returns
987 -------
988 result : `pipeBase.Struct`
989 Result of the peek.
990 Struct containing:
991 - mode : `str`
992 Peek mode that was run.
993 - binSize : `int`
994 Binning factor used.
995 - binnedSourceCat : `lsst.afw.table.SourceCatalog`
996 Source catalog from the binned exposure.
997 """
998 if mode == "auto":
999 # Note, no attempt to handle dispersed donuts. They'll default to
1000 # donut mode.
1001 if donutDiameter > self.config.donutThreshold:
1002 mode = "donut"
1003 else:
1004 if exposure.getInfo().getVisitInfo().instrumentLabel == "LATISS":
1005 # only LATISS images *can* be dispersed, and isDispersedExp
1006 # only works cleanly for LATISS
1007 mode = "spec" if isDispersedExp(exposure) else "photo"
1008 else:
1009 mode = "photo"
1011 match mode:
1012 case "donut":
1013 result = self.donut.run(exposure, donutDiameter, binSize=binSize)
1014 binSizeOut = result.binSize
1015 case "spec":
1016 result = self.spec.run(exposure, binSize=binSize)
1017 binSizeOut = result.binSize
1018 if len(result.binnedSourceCat) == 0:
1019 self.log.warning("No sources found in spec mode.")
1020 if self.config.doPhotoFallback:
1021 self.log.warning("Falling back to photo mode.")
1022 # Note that spec.run already rebinned the image,
1023 # so we don't need to do it again.
1024 result = self.photo.run(exposure, binSize=1)
1025 case "photo":
1026 result = self.photo.run(exposure, binSize=binSize)
1027 binSizeOut = result.binSize
1028 case _:
1029 raise ValueError(f"Unknown mode {mode}")
1030 return result.mode, binSizeOut, result.binnedSourceCat
1032 def transformTable(self, binSize: int, binnedSourceCat: afwTable.SourceCatalog) -> astropy.table.Table:
1033 """Make an astropy table from the source catalog but with
1034 transformations back to the original unbinned coordinates.
1036 Since there's some ambiguity in the apFlux apertures when binning,
1037 we'll only populate the table with the slots columns (slot_apFlux
1038 doesn't indicate an aperture radius). For simplicity, do the same for
1039 centroids and shapes too.
1041 And since we're only copying over the slots_* columns, we remove the
1042 "slots_" part of the column names and lowercase the first remaining
1043 letter.
1045 Parameters
1046 ----------
1047 binSize : `int`
1048 Binning factor used.
1049 binnedSourceCat : `lsst.afw.table.SourceCatalog`
1050 Source catalog from the binned exposure.
1052 Returns
1053 -------
1054 table : `astropy.table.Table`
1055 Curated source table in unbinned coordinates.
1056 """
1057 table = binnedSourceCat.asAstropy()
1058 cols = [n for n in table.colnames if n.startswith("slot")]
1059 table = table[cols]
1060 if "slot_Centroid_x" in cols:
1061 table["slot_Centroid_x"] = binSize * table["slot_Centroid_x"] + (binSize - 1) / 2
1062 table["slot_Centroid_y"] = binSize * table["slot_Centroid_y"] + (binSize - 1) / 2
1063 if "slot_Shape_x" in cols:
1064 table["slot_Shape_x"] = binSize * table["slot_Shape_x"] + (binSize - 1) / 2
1065 table["slot_Shape_y"] = binSize * table["slot_Shape_y"] + (binSize - 1) / 2
1066 table["slot_Shape_xx"] *= binSize**2
1067 table["slot_Shape_xy"] *= binSize**2
1068 table["slot_Shape_yy"] *= binSize**2
1069 # area and npixels are just confusing when binning, so remove.
1070 if "slot_PsfFlux_area" in cols:
1071 del table["slot_PsfFlux_area"]
1072 if "slot_PsfFlux_npixels" in cols:
1073 del table["slot_PsfFlux_npixels"]
1075 table.rename_columns(
1076 [n for n in table.colnames if n.startswith("slot_")],
1077 [n[5:6].lower() + n[6:] for n in table.colnames if n.startswith("slot_")],
1078 )
1080 return table
1082 def getBrightest(
1083 self, binnedSourceCat: afwTable.SourceCatalog, binSize: int, goodSourceMask: npt.NDArray[np.bool_]
1084 ) -> tuple[int, geom.Point2D, afwGeom.Quadrupole]:
1085 """Find the brightest source in the catalog.
1087 Parameters
1088 ----------
1089 binnedSourceCat : `lsst.afw.table.SourceCatalog`
1090 Source catalog from the binned exposure.
1091 binSize : `int`
1092 Binning factor used.
1093 goodSourceMask : `numpy.ndarray`
1094 Boolean array indicating which sources are good.
1096 Returns
1097 -------
1098 maxFluxIdx : `int`
1099 Index of the brightest source in the catalog.
1100 brightCentroid : `lsst.geom.Point2D`
1101 Centroid of the brightest source (unbinned coords).
1102 brightShape : `lsst.afw.geom.Quadrupole`
1103 Shape of the brightest source (unbinned coords).
1104 """
1105 fluxes = np.array([source.getApInstFlux() for source in binnedSourceCat])
1106 idxs = np.arange(len(binnedSourceCat))
1108 good = goodSourceMask & np.isfinite(fluxes)
1110 if np.sum(good) == 0:
1111 maxFluxIdx = IDX_SENTINEL
1112 brightCentroid = Point2D(np.nan, np.nan)
1113 brightShape = Quadrupole(np.nan, np.nan, np.nan)
1114 return maxFluxIdx, brightCentroid, brightShape
1116 fluxes = fluxes[good]
1117 idxs = idxs[good]
1118 maxFluxIdx = idxs[np.nanargmax(fluxes)]
1119 brightest = binnedSourceCat[maxFluxIdx]
1121 # Convert binned coordinates back to original unbinned
1122 # coordinates
1123 brightX, brightY = brightest.getCentroid()
1124 brightX = binSize * brightX + (binSize - 1) / 2
1125 brightY = binSize * brightY + (binSize - 1) / 2
1126 brightCentroid = Point2D(brightX, brightY)
1127 brightIXX = brightest.getIxx() * binSize**2
1128 brightIXY = brightest.getIxy() * binSize**2
1129 brightIYY = brightest.getIyy() * binSize**2
1130 brightShape = Quadrupole(brightIXX, brightIYY, brightIXY)
1132 return maxFluxIdx, brightCentroid, brightShape
1134 def getPsfShape(
1135 self, binnedSourceCat: afwTable.SourceCatalog, binSize: int, goodSourceMask: npt.NDArray[np.bool_]
1136 ) -> afwGeom.Quadrupole:
1137 """Estimate the modal PSF shape from the sources.
1139 Parameters
1140 ----------
1141 binnedSourceCat : `lsst.afw.table.SourceCatalog`
1142 Source catalog from the binned exposure.
1143 binSize : `int`
1144 Binning factor used.
1145 goodSourceMask : `numpy.ndarray`
1146 Boolean array indicating which sources are good.
1148 Returns
1149 -------
1150 psfShape : `lsst.afw.geom.Quadrupole`
1151 Estimated PSF shape (unbinned coords).
1152 """
1153 fluxes = np.array([source.getApInstFlux() for source in binnedSourceCat])
1154 idxs = np.arange(len(binnedSourceCat))
1156 good = goodSourceMask & np.isfinite(fluxes)
1158 if np.sum(good) == 0:
1159 return Quadrupole(np.nan, np.nan, np.nan)
1161 fluxes = fluxes[good]
1162 idxs = idxs[good]
1164 psfIXX = _estimateMode(np.array([source.getIxx() for source in binnedSourceCat])[goodSourceMask])
1165 psfIYY = _estimateMode(np.array([source.getIyy() for source in binnedSourceCat])[goodSourceMask])
1166 psfIXY = _estimateMode(np.array([source.getIxy() for source in binnedSourceCat])[goodSourceMask])
1168 return Quadrupole(
1169 psfIXX * binSize**2,
1170 psfIYY * binSize**2,
1171 psfIXY * binSize**2,
1172 )
1174 def transformShapes(
1175 self, shapes: afwGeom.Quadrupole, exposure: afwImage.Exposure, binSize: int
1176 ) -> tuple[list[afwGeom.Quadrupole], list[afwGeom.Quadrupole]]:
1177 """Transform shapes from x/y pixel coordinates to equitorial and
1178 horizon coordinates.
1180 Parameters
1181 ----------
1182 shapes : `list` of `lsst.afw.geom.Quadrupole`
1183 List of shapes (in pixel coordinates) to transform.
1184 exposure : `lsst.afw.image.Exposure`
1185 Exposure containing WCS and VisitInfo for transformation.
1186 binSize : `int`
1187 Binning factor used.
1189 Returns
1190 -------
1191 equatorialShapes : `list` of `lsst.afw.geom.Quadrupole`
1192 List of shapes transformed to equitorial (North and West)
1193 coordinates. Units are arcseconds.
1194 altAzShapes : `list` of `lsst.afw.geom.Quadrupole`
1195 List of shapes transformed to alt/az coordinates. Units are
1196 arcseconds.
1197 """
1198 pt = Point2D(np.array([*exposure.getBBox().getCenter()]) / binSize)
1199 wcs = exposure.wcs
1200 visitInfo = exposure.info.getVisitInfo()
1201 parAngle = visitInfo.boresightParAngle
1203 equatorialShapes = []
1204 altAzShapes = []
1205 for shape in shapes:
1206 if wcs is None:
1207 equatorialShapes.append(Quadrupole(np.nan, np.nan, np.nan))
1208 altAzShapes.append(Quadrupole(np.nan, np.nan, np.nan))
1209 continue
1210 # The WCS transforms to N (dec) and E (ra), but we want N and W to
1211 # conform with weak-lensing conventions. So we flip the [0]
1212 # component of the transformation.
1213 neTransform = wcs.linearizePixelToSky(pt, arcseconds).getLinear()
1214 nwTransform = LinearTransform(np.array([[-1, 0], [0, 1]]) @ neTransform.getMatrix())
1215 equatorialShapes.append(shape.transform(nwTransform))
1217 # To get from N/W to alt/az, we need to additionally rotate by the
1218 # parallactic angle.
1219 rot = LinearTransform.makeRotation(parAngle).getMatrix()
1220 aaTransform = LinearTransform(nwTransform.getMatrix() @ rot)
1221 altAzShapes.append(shape.transform(aaTransform))
1223 return equatorialShapes, altAzShapes
1225 def updateDisplay(
1226 self,
1227 exposure: afwImage.Exposure,
1228 binSize: int,
1229 binnedSourceCat: afwTable.SourceCatalog,
1230 maxFluxIdx: int,
1231 doDisplayIndices: bool,
1232 ) -> None:
1233 """Update the afwDisplay with the exposure and sources.
1235 Parameters
1236 ----------
1237 exposure : `lsst.afw.image.Exposure`
1238 Exposure to peek.
1239 binSize : `int`
1240 Binning factor used.
1241 binnedSourceCat : `lsst.afw.table.SourceCatalog`
1242 Source catalog from the binned exposure.
1243 maxFluxIdx : `int`
1244 Index of the brightest source in the catalog.
1245 doDisplayIndices : `bool`
1246 Display the source indices?
1247 """
1248 if self._display is None:
1249 raise RuntimeError("Display failed as no display provided during init()")
1251 visitInfo = exposure.info.getVisitInfo()
1252 self._display.mtv(exposure)
1253 wcs = exposure.wcs
1254 if wcs is not None:
1255 plotRose(
1256 self._display,
1257 wcs,
1258 Point2D(200 / binSize, 200 / binSize),
1259 parAng=visitInfo.boresightParAngle,
1260 len=100 / binSize,
1261 )
1263 for idx, source in enumerate(binnedSourceCat):
1264 x, y = source.getCentroid()
1265 sh = source.getShape()
1266 self._display.dot(sh, x, y)
1267 if doDisplayIndices:
1268 self._display.dot(str(idx), x, y)
1270 if maxFluxIdx != IDX_SENTINEL:
1271 self._display.dot(
1272 "+",
1273 *binnedSourceCat[maxFluxIdx].getCentroid(),
1274 ctype=afwDisplay.RED,
1275 size=10,
1276 )