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