Coverage for python/lsst/summit/utils/plotRadialAnalysis.py: 0%
239 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-09-08 02:26 -0700
1# This file is part of summit_utils.
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/>.
22from __future__ import annotations
24from collections.abc import Iterable
25from typing import Any
27import matplotlib
28import numpy as np
29import pandas
30from scipy.optimize import curve_fit # type: ignore
32from lsst.afw.cameraGeom import FIELD_ANGLE, Detector
33from lsst.daf.butler import Butler, DatasetRef
34from lsst.geom import Box2I, Extent2I, Point2I
35from lsst.summit.utils.utils import getCameraFromInstrumentName
36from lsst.utils.plotting.figures import make_figure
39def gaussian2dFitFunction(
40 xy: tuple[np.ndarray, np.ndarray],
41 peak: float,
42 fwhm: float,
43 x0: float,
44 y0: float,
45 baseline: float = 0.0,
46) -> np.ndarray:
47 """Gaussian distribution with centroid.
49 Parameters
50 ----------
51 xy: `tuple` of `np.ndarray`
52 Points coordinates.
53 peak: `float`
54 Values of the intesity peak.
55 fwhm: `float`
56 Full Width at Half Maximum fo the distribution.
57 x0: `float`
58 The x position of the 2d Guassian function.
59 y0: `float`
60 The y position of the 2d Guassian function.
61 baseline: `float`
62 Offset to apply. Default 0.
64 Returns
65 -------
66 pdf: `np.ndarray`
67 Probability density function of the distribution.
68 """
69 x, y = xy
70 sigma = fwhm / (2 * np.sqrt(2 * np.log(2)))
71 return baseline + peak * np.exp(-((x - x0) ** 2 + (y - y0) ** 2) / (2 * sigma**2))
74def moffat2dFitFunction(
75 xy: tuple[np.ndarray, np.ndarray],
76 peak: float,
77 alpha: float,
78 beta: float,
79 x0: float,
80 y0: float,
81 baseline: float,
82) -> np.ndarray:
83 """Moffat distribution with centroid.
85 Parameters
86 ----------
87 xy: `tuple` of `np.ndarray`
88 Points coordinates.
89 peak: `float`
90 Values of the intesity peak.
91 alpha: `float`
92 The alpha parameter for the Moffat distribution.
93 beta: `float`
94 The beta parameter for the Moffat distribution.
95 x0: `float`
96 x coordinate of the distribution.
97 y0: `float`
98 y coordinate of the distribution.
99 baseline: `float`
100 Offset to apply. Default 0.
102 Returns
103 -------
104 pdf: `np.ndarray`
105 Probability density function of the distribution.
106 """
107 x, y = xy
108 return baseline + peak * (1 + (((x - x0) ** 2 + (y - y0) ** 2)) / alpha**2) ** (-beta)
111def doRadialAnalysis(
112 data: np.ndarray, fitModel: str
113) -> tuple[np.ndarray, np.ndarray, np.ndarray, float, float, float]:
114 """Perform the radial analysis on a star cutout
116 Parameters
117 ----------
118 data: `np.ndarray`
119 2D array containing the star cutout.
120 fitModel: `str`
121 Model used for the fit ('moffat' or 'gauss').
123 Returns
124 -------
125 x: `np.ndarray`
126 1d array with the radial from the centroid.
127 yScatter: `np.ndarray`
128 The flattened radial distribution of the start intensity.
129 Usefull for plotting purpose.
130 yFit: `np.ndarray`
131 The fitted radial distribution.
132 fwhmFit: `float`
133 The Full Width at Half Maximum of the fitted distribution.
134 eE50Diameter: `float`
135 The Encircled Energy diameter at 50%.
136 eE80Diameter: `float`
137 The Encircled Energy diameter at 80%.
138 """
139 # Create meshgrid for fitting with x and y positions
140 xGrid, yGrid = np.meshgrid(np.arange(data.shape[1]), np.arange(data.shape[0]))
141 xy = (xGrid.ravel(), yGrid.ravel())
142 radialValues = data.ravel()
144 match fitModel.lower():
145 case "moffat":
146 # Initial guess for Moffat fitting
147 initialGuess = [
148 np.max(radialValues),
149 4,
150 2,
151 data.shape[1] / 2,
152 data.shape[0] / 2,
153 np.median(radialValues),
154 ]
155 params, _ = curve_fit(moffat2dFitFunction, xy, radialValues, p0=initialGuess, maxfev=10000)
156 _, alphaFit, betaFit, x0Fit, y0Fit, baselineFit = params
157 fwhmFit = np.abs(2.0 * alphaFit * np.sqrt(2 ** (1 / betaFit) - 1.0))
158 case "gauss":
159 # Initial guess for Gaussian fitting
160 initialGuess = [
161 np.max(radialValues),
162 10,
163 data.shape[1] / 2,
164 data.shape[0] / 2,
165 np.median(radialValues),
166 ]
167 params, _ = curve_fit(gaussian2dFitFunction, xy, radialValues, p0=initialGuess, maxfev=10000)
168 _, fwhmFit, x0Fit, y0Fit, baselineFit = params
169 case _:
170 raise ValueError(f"The model {fitModel} is not among the available ones (gauss, moffat)")
172 # Compute the curve of growth (cumulative energy)
173 radii = np.sqrt((xGrid - x0Fit) ** 2 + (yGrid - y0Fit) ** 2).ravel()
175 sortedIndices = np.argsort(radii)
176 sortedradii = radii[sortedIndices]
177 sortedValues = radialValues[sortedIndices] - baselineFit
178 cumulativeEnergy = np.cumsum(sortedValues)
180 # Determine 50% and 80% encircled energy diameters
181 eE50Diameter = 2 * sortedradii[np.searchsorted(cumulativeEnergy, 0.5 * cumulativeEnergy[-1])]
182 eE80Diameter = 2 * sortedradii[np.searchsorted(cumulativeEnergy, 0.8 * cumulativeEnergy[-1])]
184 x = 0.2 * sortedradii
185 yScatter = sortedValues + baselineFit # added back the backgroud.
187 match fitModel.lower():
188 case "moffat":
189 yFit = moffat2dFitFunction(xy, *params)[sortedIndices]
190 case "gauss":
191 yFit = gaussian2dFitFunction(xy, *params)[sortedIndices]
192 case _:
193 raise ValueError(f"The model {fitModel} is not among the available ones (gauss, moffat)")
195 return (
196 x,
197 yScatter,
198 yFit,
199 fwhmFit,
200 eE50Diameter,
201 eE80Diameter,
202 )
205def makeLayerPlot(
206 ax: matplotlib.axes.Axes,
207 data: np.ndarray,
208 e1: float,
209 e2: float,
210 e: float,
211 fitModel: str,
212 layers: list[str] | str = "all",
213 levels: np.ndarray | Iterable[float] | None = None,
214) -> tuple[float, float, float]:
215 """Make per axes layer plot.
217 Create a plot with three possible layers:
218 - Background image with the star cutout
219 - The contour level of the star
220 - The radial profile with Gaussian and Moffat fit
221 The value of FWHM and Encircled Energy Radii (EE)
222 at 50% and 80% are reported if the 'radial' layer has been chosen.
224 Parameters
225 ----------
226 ax: `matplotlib.axes.Axes`
227 The axes to use
228 data: `np.ndarray`
229 2D array containing the star cutout
230 fitModel: `str`
231 Model used for the fit ('moffat' or 'gauss')
232 layers: list[str] | str = "all",
233 List of layers to be displayed
234 ('background', 'radial', 'contour', 'ellipticity').
235 levels: `np.ndarray` or `Iterable` of `float` or `None`, optional
236 The levels value for the contour layer.
237 If None, is set to `np.linspace(1.5*np.std(data), data.max(), 5)`
239 Returns
240 -------
241 fwhmFit: `float`
242 The FWHM of the fitted distribution.
243 eE50Diameter: `float`
244 The encircled energy diameter at 50%.
245 eE80Diameter: `float`
246 The encircled energy diameter at 80%.
247 """
248 if layers == "all":
249 layers = ["background", "radial", "contours", "ellipticity"]
251 (
252 x,
253 yScatter,
254 yFit,
255 fwhmFit,
256 eE50Diameter,
257 eE80Diameter,
258 ) = doRadialAnalysis(data, fitModel)
260 # get figure and axes position on figure
261 # to create multiple axes on that position
262 fig = ax.get_figure()
263 bbox = ax.get_position()
264 bboxArray: tuple[float, float, float, float] = bbox.bounds
266 assert fig is not None
268 # plot the background layer if present
269 axBkg = None
270 if "background" in layers:
271 axBkg = fig.add_axes(bboxArray, frameon=False)
272 axBkg.imshow(data, cmap="gray", origin="lower", zorder=1)
273 axBkg.set(zorder=1, xticks=[], yticks=[])
275 # plot the contour layer if present
276 if "contour" in layers:
277 axCtr = fig.add_axes(bboxArray, frameon=False)
278 xGrid, yGrid = np.meshgrid(np.arange(data.shape[1]), np.arange(data.shape[0]))
279 levels = levels if levels is not None else np.linspace(1.5 * np.std(data), data.max(), 5)
280 axCtr.contour(xGrid, yGrid, data, cmap="spring", levels=levels, alpha=0.7)
281 axCtr.set(zorder=2, facecolor=(1, 1, 1, 0), xticks=[], yticks=[])
283 # plot the radial analysis layer if present
284 if "radial" in layers:
285 axRad = fig.add_axes(bboxArray, frameon=False)
286 axRad.scatter(
287 x,
288 yScatter,
289 marker="o",
290 s=5,
291 facecolor="none",
292 edgecolor="cyan",
293 label="Data",
294 )
296 axRad.plot(
297 x,
298 yFit,
299 ls="-",
300 color="y",
301 label=f"{fitModel} Fit",
302 )
304 axRad.set(zorder=3, facecolor=(1, 1, 1, 0), xticks=[], yticks=[])
306 if "ellipticity" in layers:
307 axEll = fig.add_axes(bboxArray, frameon=False)
308 axEll.set(zorder=4, facecolor=(1, 1, 1, 0), xticks=[], yticks=[])
310 quiverScale = 1
311 eU = e * np.cos(0.5 * np.arctan2(e2, e1))
312 eV = e * np.sin(0.5 * np.arctan2(e2, e1))
313 peak = np.argmax(data)
314 centerY, centerX = np.unravel_index(peak, data.shape)
316 # better to use the background axes limits
317 if "background" in layers:
318 assert axBkg is not None, "Logic error: axBkg should have been defined above"
319 axEll.set(xlim=axBkg.get_xlim(), ylim=axBkg.get_ylim())
320 else: # otherwise use the cutout limits
321 axEll.set(xlim=(0, data.shape[1]), ylim=(0, data.shape[0]))
323 axEll.quiver(
324 centerX,
325 centerY,
326 eU,
327 eV,
328 headlength=0,
329 headaxislength=0,
330 scale=quiverScale,
331 pivot="mid",
332 color="fuchsia",
333 width=0.015,
334 )
336 return fwhmFit, eE50Diameter, eE80Diameter
339def compactifyLayout(
340 rectDict: dict[str, tuple[float, float, float, float]],
341) -> dict[str, tuple[float, float, float, float]]:
342 """Compact the layout of the rectangles to fit in a smaller area.
343 This function rescales the rectangles to fit within a specified area
344 while maintaining their aspect ratios.
346 Parameters
347 ----------
348 rectDict: `dict`
349 Dictionary with the detector boxes.
351 Returns
352 -------
353 rectDict: `dict`
354 Dictionary with the rescaled rectangles.
355 """
357 margin = 0.02
359 lefts = [r[0] for r in rectDict.values()]
360 bottoms = [r[1] for r in rectDict.values()]
362 min_left = min(lefts)
363 max_left = max(lefts)
364 min_bottom = min(bottoms)
365 max_bottom = max(bottoms)
367 range_x = max_left - min_left if max_left != min_left else 1
368 range_y = max_bottom - min_bottom if max_bottom != min_bottom else 1
370 scale_x = (1 - 2 * margin) / range_x
371 scale_y = (1 - 2 * margin) / range_y
373 newRectDict: dict[str, tuple[float, float, float, float]] = {}
374 for name, (left, bottom, width, height) in rectDict.items():
375 new_left = (left - min_left) * scale_x + margin
376 new_bottom = (bottom - min_bottom) * scale_y + margin
377 newRectDict[name] = (new_left, new_bottom, width, height)
379 return newRectDict
382def computeColorbarRect(
383 rectDict: dict[str, tuple[float, float, float, float]],
384) -> tuple[float, float, float, float]:
385 """Compute the colorbar rectangle based on the detector rectangles
387 Parameters
388 ----------
389 rectDict: `dict`
390 Dictionary with the detector boxes.
392 Returns
393 -------
394 rect: `tuple`
395 The box in which the colorbar will be drawn.
396 """
398 width = 0.015
399 padding = 0.01
401 rects = list(rectDict.values())
403 min_bottom = min(r[1] for r in rects)
404 max_top = max(r[1] + r[3] for r in rects)
405 max_right = max(r[0] + r[2] for r in rects)
407 height = max_top - min_bottom
408 left = max_right + padding
409 bottom = min_bottom
411 return (left, bottom, width, height)
414def createFigWithInstrumentLayout(
415 fig: matplotlib.figure.Figure,
416 instrument: str,
417 onlyS11: bool = False,
418) -> dict[str, matplotlib.axes.Axes]:
419 """Create a figure with the requested instrument layout
421 Parameters
422 ----------
423 fig: `matplotlib.figure.Figure`
424 The figure to use
425 instrument: `str`
426 The instrument name.
427 onlyS11: `bool`, optional
428 If True, only S11 detectors are shown. Default False.
430 Returns
431 -------
432 axsDict: `dict[str, matplotlib.axes.Axes]`
433 A dictionary with the detector names as keys and the axes as values.
434 """
436 camera = getCameraFromInstrumentName(instrument)
437 detectors = [detector.getId() for detector in camera]
439 rectDict: dict[str, tuple[float, float, float, float]] = {}
440 for name in detectors:
441 detector = camera.get(name)
442 detName = detector.getName()
444 # Get the corners of the detector in FIELD_ANGLE
445 corners = detector.getCorners(FIELD_ANGLE)
446 corners_deg = np.rad2deg(corners) # Convert corners to degrees
448 # turn the cornsers coords in (letf, bottom, width, height)
449 detRect = (
450 corners_deg[:, 0].min(),
451 corners_deg[:, 1].min(),
452 corners_deg[:, 0].max() - corners_deg[:, 0].min(),
453 corners_deg[:, 1].max() - corners_deg[:, 1].min(),
454 )
456 if onlyS11 and "S11" not in detName:
457 continue
458 rectDict[detName] = detRect
460 if onlyS11:
461 rectDict = compactifyLayout(rectDict)
463 axsDict = {}
464 for detName, rect in rectDict.items():
465 # Create an axes for each detector
466 ax = fig.add_axes(rect)
467 ax.set(xticks=[], yticks=[])
468 axsDict[detName] = ax
470 cbraRect = computeColorbarRect(rectDict)
471 cbarAx = fig.add_axes(cbraRect)
472 axsDict["cbar"] = cbarAx
474 return axsDict
477def makePsfPanel(
478 cutouts: dict[str, tuple[np.ndarray, tuple[float, float, float]]],
479 instrument: str = "LSSTComCam",
480 onlyS11: bool = False,
481 layers: list[str] | str = "all",
482 fitModel: str = "moffat",
483 levels: np.ndarray | Iterable[float] | None = None,
484 **kwargs: Any,
485) -> matplotlib.figure.Figure:
486 """Make a per-detector PSF radial analysis.
488 Each subplot shows for a detector a PSF cutout, a radial analysis and the
489 morphology contour, a custom selection of this layer is possible.
491 Parameters
492 ----------
493 cutouts: `dict[str, np.ndarray]`
494 A detector's name key dictionary containing
495 the 2D array of the star cutouts.
496 instrument: `str`, optional
497 Detector type. Default 'LSSTComCam'.
498 onlyS11: `bool`, optional
499 If True, only S11 detectors are shown. Default False.
500 layers: `str` or `list` of `str`, optional
501 List of layers to be displayed ('background', 'radial', 'contour').
502 It is possible to pass also the string 'all' as a shortcut for
503 ['background', 'radial', 'contour']. Default 'all'.
504 fitModel: `str`, optional
505 Model used for the radial fit ('moffat' or 'gauss').
506 Default 'moffat'.
507 levels: `np.ndarray` or `Iterable` of `float` or `None`, optional
508 The levels value for the contour layer.
509 If None, is set to `np.linspace(1.5*np.std(data), data.max(), 5)`.
510 **kwargs:
511 Additional keyword arguments passed to matplotlib Figure constructor.
513 Returns
514 -------
515 fig: `matplotlib.figure.Figure`
516 The figure.
517 """
518 fig = make_figure(**kwargs)
519 axsDict = createFigWithInstrumentLayout(fig, instrument, onlyS11=onlyS11)
521 if layers == "all":
522 layers = ["background", "radial", "contours", "ellipticity"]
524 fwhmDict: dict[str, float] = {}
525 ee50Dict: dict[str, float] = {}
526 ee80Dict: dict[str, float] = {}
527 for detName, (cutout, (e1, e2, e)) in cutouts.items():
528 fwhm, ee50, ee80 = makeLayerPlot(axsDict[detName], cutout, e1, e2, e, fitModel, layers, levels)
529 fwhmDict[detName] = fwhm
530 ee50Dict[detName] = ee50
531 ee80Dict[detName] = ee80
533 cmap = matplotlib.colormaps["RdYlGn_r"]
534 for detName, (_, (e1, e2, _)) in cutouts.items():
535 bbox = axsDict[detName].get_position()
536 bboxArray = bbox.bounds
538 axText = fig.add_axes(bboxArray, frameon=True)
539 axText.set(zorder=4, facecolor=(1, 1, 1, 0), xticks=[], yticks=[])
541 if detName == "cbar":
542 continue
544 val = (fwhmDict[detName] - min(fwhmDict.values())) / (max(fwhmDict.values()) - min(fwhmDict.values()))
545 color = cmap(val)
546 axText.patch.set(lw=7, ec=color)
547 for _, spine in axText.spines.items():
548 spine.set_color(color)
550 # Text with FWHM and EE50|80
551 text = (
552 f'FWHM: {fwhmDict[detName] * 0.2:.2f}" '
553 f'EE80: {ee80Dict[detName] * 0.2:.2f}"\n'
554 f"E1|2: {e1:.2f}|{e2:.2f}"
555 )
556 axText.text(
557 0.97,
558 0.95,
559 text,
560 color=color,
561 fontsize=12,
562 fontweight="bold",
563 transform=axText.transAxes,
564 ha="right",
565 va="top",
566 )
568 axText.text(
569 0.3,
570 0.1,
571 detName,
572 color="silver",
573 fontsize=11,
574 fontweight="bold",
575 transform=axText.transAxes,
576 ha="right",
577 va="top",
578 )
580 # set colorbar
581 vmin = None
582 vmax = None
583 if fwhmDict:
584 vmin = min(fwhmDict.values()) * 0.2
585 vmax = max(fwhmDict.values()) * 0.2
586 cbar = fig.colorbar(
587 matplotlib.cm.ScalarMappable(
588 cmap=cmap,
589 norm=matplotlib.pyplot.Normalize(vmin=vmin, vmax=vmax),
590 ),
591 cax=axsDict["cbar"],
592 )
593 cbar.ax.set_ylabel(ylabel='FWHM"', fontsize=30)
594 cbar.ax.tick_params(labelsize=30)
596 return fig
599def generateCutout(
600 butler: Butler,
601 imgRef: DatasetRef,
602 detector: Detector,
603 target: np.ndarray | list[float] | tuple[float, float],
604) -> np.ndarray:
605 """Generate the cutout around a target position
607 Parameters
608 ----------
609 butler: `lsst.daf.butler.Butler`
610 The butler to use to get the image
611 imgRef: `lsst.daf.butler.DatasetRef`
612 The dataset reference to use to get the image
613 detector: `lsst.afw.cameraGeom.Detector`
614 The detector to use to get calculate the region of interest
615 target: `np.ndarray` or `list` of `float` or `tuple` of `float`
616 The coordinates of the cutout center
618 Returns
619 -------
620 cutout: `np.ndarray`
621 The square cutout around the center position.
622 """
624 pad = 20
625 detBbox = detector.getBBox()
626 start = Point2I(target[0] - (pad // 2), target[1] - (pad // 2))
627 dim = Extent2I(pad, pad)
628 roiBbox = detBbox.clippedTo(Box2I(start, dim))
629 cutout = butler.get(imgRef, parameters={"bbox": roiBbox}).image.array
631 return cutout
634def findNearestStarToCenter(
635 tab: pandas.DataFrame,
636 detector: Detector,
637 instrument: str,
638) -> tuple[np.ndarray, tuple[float, float, float]]:
639 """Find the nearest star w.r.t to the detector center
640 N.B. The seacrh is done in PIXEL coordinates.
642 Parameters
643 ----------
644 tab: `pandas.DataFrame`
645 pandas.DataFrame with the in focus stars positions.
646 detector: `lsst.afw.cameraGeom.Detector`
647 The detector realted to the sourceTable's positions.
648 instrument: `str`
649 Instrument name.
651 Returns
652 -------
653 `np.ndarray`
654 The x and y coordinates of the nearest star.
655 `tuple`
656 The ellipticity parameters (e1, e2, e).
657 """
659 if instrument == "LSSTComCam":
660 xCol = "slot_Centroid_x"
661 yCol = "slot_Centroid_y"
662 else: # for now just work with src file that has the same column.
663 xCol = "slot_Centroid_x" # "x"
664 yCol = "slot_Centroid_y" # "y"
666 target = (detector.getBBox().centerX, detector.getBBox().centerY)
668 tab["center_sep"] = np.sqrt((tab[xCol] - target[0]) ** 2 + (tab[yCol] - target[1]) ** 2)
669 most_close = tab.sort_values(by=["center_sep"]).iloc[0].name
670 nearest = tab.loc[most_close, [xCol, yCol]].values
672 # from makeTableFromSourceCatalogs on lsst.summit_extras.plotting
673 iXX = tab.loc[most_close, "slot_Shape_xx"] * (0.2) ** 2
674 iYY = tab.loc[most_close, "slot_Shape_yy"] * (0.2) ** 2
675 iXY = tab.loc[most_close, "slot_Shape_xy"] * (0.2) ** 2
676 T = iXX + iYY
677 e1 = (iXX - iYY) / T
678 e2 = 2 * iXY / T
679 e = np.hypot(e1, e2)
681 return nearest, (e1, e2, e)
684def makePanel(
685 butler: Butler,
686 visit: int,
687 onlyS11: bool = False,
688 **kwargs: Any,
689) -> matplotlib.figure.Figure:
690 """Create the panel with the in focus stars.
691 See the documentation of `makePsfPanel` for more information.
693 Parameters
694 ----------
695 butler: `lsst.daf.butler.Butler`
696 The butler to use to get the image and the source table
697 imgRefs: `lsst.daf.butler.DatasetRef`
698 The dataset reference to use to get the image
699 srcRefs: `lsst.daf.butler.DatasetRef`
700 The dataset reference to use to get the source table
701 onlyS11: `bool`, optional
702 If True, only S11 detectors are shown. Default False.
703 **kwargs:
704 Parameters for the `makePsfPanel` method.
706 Returns
707 -------
708 fig: `matplotlib.figure.Figure`
709 The figure.
711 Raises
712 ------
713 ValueError
714 If no image or source table datasets are found for the given visit.
715 """
717 # retrieve the image and source table dataset references
718 imgRefs = butler.query_datasets("post_isr_image", where=f"exposure={visit}", explain=False)
719 srcRefs = butler.query_datasets("single_visit_psf_star", where=f"exposure={visit}", explain=False)
720 if not imgRefs or not srcRefs:
721 raise ValueError(f"No image and source tables found for visit {visit}")
723 # grab the instrument name from one of the imgRefs
724 instrument = imgRefs[0].dataId["instrument"]
725 assert isinstance(instrument, str), f"Instrument name {instrument} is not a string"
726 camera = getCameraFromInstrumentName(instrument)
727 instrumentName = camera.getName()
729 # if only S11 then keep only the S11 detectors
730 if onlyS11:
731 imgRefs = [dr for dr in imgRefs if "S11" in camera[dr.dataId["detector"]].getName()]
732 srcRefs = [dr for dr in srcRefs if "S11" in camera[dr.dataId["detector"]].getName()]
734 # retrieve the detector object for each detector
735 detNameDict = {detector.getName(): detector for detector in camera}
737 # interesct the detNum for images and tables
738 imgDetName = {camera[dr.dataId["detector"]].getName() for dr in imgRefs}
739 srcDetName = {camera[dr.dataId["detector"]].getName() for dr in srcRefs}
740 commonDetName = imgDetName.intersection(srcDetName)
742 # first retrieve the srcTable datasets from butler
743 filterColumn = "calib_psf_used"
744 sourceTableDict = {
745 camera[dr.dataId["detector"]].getName(): butler.get(dr).to_pandas()
746 for dr in srcRefs
747 if camera[dr.dataId["detector"]].getName() in commonDetName
748 }
749 sourceTableDict = {detName: tab[tab[filterColumn]] for detName, tab in sourceTableDict.items()}
751 # filter commoDetName to keep only srcTable with non zero rows
752 filterDetName = []
753 for detName in commonDetName:
754 if sourceTableDict[detName].shape[0] > 0:
755 filterDetName.append(detName)
757 # find the most center star in the srcTables
758 candidates = {
759 detName: findNearestStarToCenter(sourceTableDict[detName], detNameDict[detName], instrumentName)
760 for detName in filterDetName
761 }
763 # filter the imgRefs
764 filterImgRefDict = {}
765 for imgRef in imgRefs:
766 detName = camera[imgRef.dataId["detector"]].getName()
767 if detName in filterDetName:
768 filterImgRefDict[detName] = imgRef
770 # generate the cutouts
771 cutouts = {
772 detName: (
773 generateCutout(butler, filterImgRefDict[detName], detNameDict[detName], candidates[detName][0]),
774 candidates[detName][1],
775 )
776 for detName in filterDetName
777 }
779 fig = makePsfPanel(cutouts, instrumentName, onlyS11=onlyS11, **kwargs)
780 return fig