Coverage for python / lsst / analysis / tools / actions / plot / focalPlanePlot.py: 13%
244 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-17 09:36 +0000
1# This file is part of analysis_tools.
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
24__all__ = ("FocalPlanePlot", "FocalPlaneGeometryPlot")
26from collections.abc import Mapping
28import matplotlib.patheffects as pathEffects
29import matplotlib.pyplot as plt
30import numpy as np
31from matplotlib.collections import PatchCollection
32from matplotlib.figure import Figure
33from matplotlib.offsetbox import AnchoredText
34from matplotlib.patches import Polygon
35from mpl_toolkits.axes_grid1 import make_axes_locatable
36from scipy.stats import binned_statistic_2d, binned_statistic_dd
38from lsst.afw.cameraGeom import FOCAL_PLANE, PIXELS, Camera
39from lsst.pex.config import ChoiceField, Field
41from ...interfaces import KeyedData, KeyedDataSchema, PlotAction, Scalar, Vector
42from ...math import nanMax, nanMedian, nanMin, nanSigmaMad
43from .plotUtils import addPlotInfo, mkColormap, sortAllArrays
46class FocalPlanePlot(PlotAction):
47 """Plots the focal plane distribution of a parameter.
49 Given the detector positions in x and y, the focal plane positions are
50 calculated using the camera model. A 2d binned statistic (default is mean)
51 is then calculated and plotted for the parameter z as a function of the
52 focal plane coordinates.
53 """
55 xAxisLabel = Field[str](doc="Label to use for the x axis.", default="x (mm)", optional=True)
56 yAxisLabel = Field[str](doc="Label to use for the y axis.", default="y (mm)", optional=True)
57 zAxisLabel = Field[str](doc="Label to use for the z axis.", optional=False)
58 nBins = Field[int](
59 doc="Number of bins to use within the effective plot ranges along the spatial directions.",
60 default=200,
61 )
62 doUseAdaptiveBinning = Field[bool](
63 doc="If set to True, the number of bins is adapted to the source"
64 " density, with lower densities using fewer bins. Under these"
65 " circumstances the nBins parameter sets the minimum number of bins.",
66 default=False,
67 )
68 statistic = Field[str](
69 doc="Operation to perform in binned_statistic_2d",
70 default="mean",
71 )
72 plotMin = Field[float](
73 doc="Minimum in z-value to display in the focal plane plot and in the histogram plot, if applicable",
74 default=None,
75 optional=True,
76 )
77 plotMax = Field[float](
78 doc="Maximum in z-value to display in the focal plane plot and in the histogram plot, if applicable",
79 default=None,
80 optional=True,
81 )
82 showStats = Field[bool](doc="Show statistics for plotted data", default=True)
83 addHistogram = Field[bool](doc="Add a histogram of all input points", default=False)
84 histBins = Field[int](doc="Number of bins to use in histogram", default=30)
86 def __call__(self, data: KeyedData, **kwargs) -> Mapping[str, Figure] | Figure:
87 self._validateInput(data, **kwargs)
88 return self.makePlot(data, **kwargs)
89 # table is a dict that needs: x, y, run, skymap, filter, tract,
91 def _validateInput(self, data: KeyedData, **kwargs) -> None:
92 """NOTE currently can only check that something is not a Scalar, not
93 check that the data is consistent with Vector
94 """
95 needed = self.getInputSchema(**kwargs)
96 if remainder := {key.format(**kwargs) for key, _ in needed} - {
97 key.format(**kwargs) for key in data.keys()
98 }:
99 raise ValueError(f"Task needs keys {remainder} but they were not found in input")
100 for name, typ in needed:
101 isScalar = issubclass((colType := type(data[name.format(**kwargs)])), Scalar)
102 if isScalar and typ != Scalar:
103 raise ValueError(f"Data keyed by {name} has type {colType} but action requires type {typ}")
105 def getInputSchema(self, **kwargs) -> KeyedDataSchema:
106 base = []
107 base.append(("x", Vector))
108 base.append(("y", Vector))
109 base.append(("z", Vector))
110 base.append(("statMask", Vector))
112 return base
114 def statsAndText(self, arr, mask=None):
115 """Calculate some stats from an array and return them
116 and some text.
117 """
118 numPoints = len(arr)
119 if mask is not None:
120 arr = arr[mask]
121 med = nanMedian(arr)
122 sigMad = nanSigmaMad(arr)
124 statsText = (
125 f"Median: {med:0.2f}\n"
126 + r"$\sigma_{MAD}$: "
127 + f"{sigMad:0.2f}\n"
128 + r"n$_{points}$: "
129 + f"{numPoints}"
130 )
132 return med, sigMad, statsText
134 def _addHistogram(self, histAx, data):
135 bins = np.linspace(
136 (self.plotMin if self.plotMin else nanMin(data.astype(float))),
137 (self.plotMax if self.plotMax else nanMax(data.astype(float))),
138 self.histBins,
139 )
140 histAx.hist(data.astype(float), bins=bins)
141 histAx.set_xlabel(self.zAxisLabel)
142 histAx.set_ylabel("Bin count")
143 underflow = np.count_nonzero(data < bins[0])
144 overflow = np.count_nonzero(data > bins[-1])
145 nonfinite = np.count_nonzero(~np.isfinite(data))
146 text = f"Underflow = {underflow}\nOverflow = {overflow}\nNon-Finite = {nonfinite}"
147 anchored_text = AnchoredText(text, loc=1, pad=0.5)
148 histAx.add_artist(anchored_text)
150 def makePlot(
151 self,
152 data: KeyedData,
153 camera: Camera,
154 plotInfo: Mapping[str, str] | None = None,
155 **kwargs,
156 ) -> Figure:
157 """Prep the catalogue and then make a focalPlanePlot of the given
158 column.
160 Uses the axisLabels config options `x` and `y` to make an image, where
161 the color corresponds to the 2d binned statistic (the mean is the
162 default) applied to the `z` column. A summary panel is shown in the
163 upper right corner of the resultant plot. The code uses the
164 selectorActions to decide which points to plot and the
165 statisticSelector actions to determine which points to use for the
166 printed statistics.
168 Parameters
169 ----------
170 data : `pandas.core.frame.DataFrame`
171 The catalog to plot the points from.
172 camera : `lsst.afw.cameraGeom.Camera`
173 The camera used to map from pixel to focal plane positions.
174 plotInfo : `dict`
175 A dictionary of information about the data being plotted with keys:
177 ``"run"``
178 The output run for the plots (`str`).
179 ``"skymap"``
180 The type of skymap used for the data (`str`).
181 ``"filter"``
182 The filter used for this data (`str`).
183 ``"tract"``
184 The tract that the data comes from (`str`).
185 ``"bands"``
186 The band(s) that the data comes from (`list` of `str`).
188 Returns
189 -------
190 fig : `matplotlib.figure.Figure`
191 The resulting figure.
192 """
193 if plotInfo is None:
194 plotInfo = {}
196 # This is a little hacky but is the easiest way to solve
197 # the bands being printed on the plot correctly.
198 if "band" in plotInfo.keys():
199 plotInfo["bands"] = [plotInfo["band"]]
201 if len(data["x"]) == 0:
202 noDataFig = Figure()
203 noDataFig.text(0.3, 0.5, "No data to plot after selectors applied")
204 noDataFig = addPlotInfo(noDataFig, plotInfo)
205 return noDataFig
207 if self.addHistogram:
208 fig, [ax, histAx] = plt.subplots(1, 2, dpi=300, figsize=(12, 6), width_ratios=[3, 2])
209 else:
210 fig = plt.figure(dpi=300)
211 ax = fig.add_subplot(111)
213 detectorIds = np.unique(data["detector"])
214 focalPlane_x = np.zeros(len(data["x"]))
215 focalPlane_y = np.zeros(len(data["y"]))
216 for detectorId in detectorIds:
217 detector = camera[detectorId]
218 map = detector.getTransform(PIXELS, FOCAL_PLANE).getMapping()
220 detectorInd = data["detector"] == detectorId
221 points = np.array([data["x"][detectorInd], data["y"][detectorInd]])
223 fp_x, fp_y = map.applyForward(points)
224 focalPlane_x[detectorInd] = fp_x
225 focalPlane_y[detectorInd] = fp_y
227 if self.doUseAdaptiveBinning:
228 # Use a course 32x32 binning to determine the mean source density
229 # in regions where there are sources.
230 binsx = np.linspace(focalPlane_x.min() - 1e-5, focalPlane_x.max() + 1e-5, 33)
231 binsy = np.linspace(focalPlane_y.min() - 1e-5, focalPlane_y.max() + 1e-5, 33)
233 binnedNumSrc = np.histogram2d(focalPlane_x, focalPlane_y, bins=[binsx, binsy])[0]
234 meanSrcDensity = np.mean(binnedNumSrc, where=binnedNumSrc > 0.0)
236 numBins = int(np.round(16.0 * np.sqrt(meanSrcDensity)))
237 numBins = max(numBins, self.nBins)
238 else:
239 numBins = self.nBins
241 # Add an arbitrary small offset to bins to ensure that the minimum does
242 # not equal the maximum.
243 binsx = np.linspace(focalPlane_x.min() - 1e-5, focalPlane_x.max() + 1e-5, numBins)
244 binsy = np.linspace(focalPlane_y.min() - 1e-5, focalPlane_y.max() + 1e-5, numBins)
246 statistic, x_edge, y_edge, binnumber = binned_statistic_2d(
247 focalPlane_x, focalPlane_y, data["z"], statistic=self.statistic, bins=[binsx, binsy]
248 )
249 binExtent = [x_edge[0], x_edge[-1], y_edge[0], y_edge[-1]]
251 if self.showStats:
252 sortedArrs = sortAllArrays([data["z"], data["x"], data["y"], data["statMask"]])
253 [colorVals, xs, ys, stat] = sortedArrs
254 statMed, statMad, statsText = self.statsAndText(colorVals, mask=stat)
255 bbox = dict(facecolor="paleturquoise", alpha=0.5, edgecolor="none")
256 ax.text(0.8, 0.91, statsText, transform=fig.transFigure, fontsize=8, bbox=bbox)
258 median = nanMedian(statistic.ravel())
259 mad = nanSigmaMad(statistic.ravel())
261 vmin = self.plotMin if (self.plotMin is not None) else (median - 2 * mad)
262 vmax = self.plotMax if (self.plotMax is not None) else (median + 2 * mad)
264 plot = ax.imshow(statistic.T, extent=binExtent, vmin=vmin, vmax=vmax, origin="lower")
266 divider = make_axes_locatable(ax)
267 cax = divider.append_axes("right", size="5%", pad=0.05)
268 fig.colorbar(plot, cax=cax, extend="both")
269 text = cax.text(
270 0.5,
271 0.5,
272 self.zAxisLabel,
273 color="k",
274 rotation="vertical",
275 transform=cax.transAxes,
276 ha="center",
277 va="center",
278 fontsize=10,
279 )
280 text.set_path_effects([pathEffects.Stroke(linewidth=3, foreground="w"), pathEffects.Normal()])
281 cax.tick_params(labelsize=7)
283 ax.set_xlabel(self.xAxisLabel)
284 ax.set_ylabel(self.yAxisLabel)
285 ax.tick_params(axis="x", labelrotation=25)
286 ax.tick_params(labelsize=7)
288 ax.set_aspect("equal")
290 if self.addHistogram:
291 self._addHistogram(histAx, data["z"])
293 plt.draw()
295 # Add useful information to the plot
296 plt.subplots_adjust(left=0.05, right=0.95)
297 fig = plt.gcf()
298 if plotInfo:
299 fig = addPlotInfo(fig, plotInfo)
301 return fig
304class FocalPlaneGeometryPlot(FocalPlanePlot):
305 """Plots the focal plane distribution of a parameter in afw camera
306 geometry units: amplifiers and detectors.
308 Given the detector positions in x and y, the focal plane positions
309 are calculated using the camera model. A 2d binned statistic
310 (default is mean) is then calculated and plotted for the parameter
311 z as a function of the camera geometry segment the input points
312 fall upon.
314 The ``xAxisLabel``, ``yAxisLabel``, ``zAxisLabel``, and
315 ``statistic`` variables are inherited from the parent class.
316 """
318 level = ChoiceField[str](
319 doc="Which geometry level should values be plotted?",
320 default="amplifier",
321 allowed={
322 "amplifier": "Plot values per readout amplifier.",
323 "detector": "Plot values per detector.",
324 },
325 )
327 def getInputSchema(self, **kwargs) -> KeyedDataSchema:
328 base = []
329 base.append(("detector", Vector))
330 if self.level == "amplifier":
331 base.append(("amplifier", Vector))
332 base.append(("z", Vector))
334 return base
336 def makePlot(
337 self,
338 data: KeyedData,
339 camera: Camera,
340 plotInfo: Mapping[str, str] | None = None,
341 **kwargs,
342 ) -> Figure:
343 """Prep the catalogue and then make a focalPlanePlot of the given
344 column.
346 Uses the axisLabels config options `x` and `y` to make an image, where
347 the color corresponds to the 2d binned statistic (the mean is the
348 default) applied to the `z` column. A summary panel is shown in the
349 upper right corner of the resultant plot. The code uses the
350 selectorActions to decide which points to plot and the
351 statisticSelector actions to determine which points to use for the
352 printed statistics.
354 Parameters
355 ----------
356 data : `pandas.core.frame.DataFrame`
357 The catalog to plot the points from. This is expected to
358 have the following columns/keys:
360 ``"detector"``
361 The integer detector id for the points.
362 ``"amplifier"``
363 The string amplifier name for the points.
364 ``"z"``
365 The numerical value that will be combined via
366 ``statistic`` to the binned value.
367 ``"x"``
368 Focal plane x position, optional.
369 ``"y"``
370 Focal plane y position, optional.
371 camera : `lsst.afw.cameraGeom.Camera`
372 The camera used to map from pixel to focal plane positions.
373 plotInfo : `dict`
374 A dictionary of information about the data being plotted with keys:
376 ``"run"``
377 The output run for the plots (`str`).
378 ``"skymap"``
379 The type of skymap used for the data (`str`).
380 ``"filter"``
381 The filter used for this data (`str`).
382 ``"tract"``
383 The tract that the data comes from (`str`).
384 ``"bands"``
385 The band(s) that the data comes from (`list` of `str`).
387 Returns
388 -------
389 fig : `matplotlib.figure.Figure`
390 The resulting figure.
391 """
393 cmap = mkColormap(["midnightBlue", "lightcyan", "darkgreen"])
394 cmap.set_bad(color="none")
396 if plotInfo is None:
397 plotInfo = {}
399 if len(data["z"]) == 0:
400 noDataFig = Figure()
401 noDataFig.text(0.3, 0.5, "No data to plot after selectors applied")
402 noDataFig = addPlotInfo(noDataFig, plotInfo)
403 return noDataFig
405 if self.addHistogram:
406 fig, [ax, histAx] = plt.subplots(1, 2, dpi=300, figsize=(12, 6), width_ratios=[3, 2])
407 else:
408 fig = plt.figure(dpi=300)
409 ax = fig.add_subplot(111)
411 detectorIds = np.unique(data["detector"])
412 focalPlane_x = np.zeros(len(data["z"]))
413 focalPlane_y = np.zeros(len(data["z"]))
415 patches = []
416 values = []
418 # Plot bounding box that will be used to set the axes below.
419 plotLimit_x = [0.0, 0.0]
420 plotLimit_y = [0.0, 0.0]
422 for detectorId in detectorIds:
423 detector = camera[detectorId]
425 # We can go stright to fp coordinates.
426 corners = [(c.getX(), c.getY()) for c in detector.getCorners(FOCAL_PLANE)]
427 corners = np.array(corners)
429 # U/V coordinates represent focal plane locations.
430 minU, minV = corners.min(axis=0)
431 maxU, maxV = corners.max(axis=0)
433 # See if the plot bounding box needs to be extended:
434 if minU < plotLimit_x[0]:
435 plotLimit_x[0] = minU
436 if minV < plotLimit_y[0]:
437 plotLimit_y[0] = minV
438 if maxU > plotLimit_x[1]:
439 plotLimit_x[1] = maxU
440 if maxV > plotLimit_y[1]:
441 plotLimit_y[1] = maxV
443 # X/Y coordinates represent detector internal coordinates.
444 # Detector extent in detector coordinates
445 minX, minY = detector.getBBox().getMin()
446 maxX, maxY = detector.getBBox().getMax()
448 if self.level.lower() == "detector":
449 detectorInd = data["detector"] == detectorId
451 # This does the appropriate statistic for this
452 # detector's data.
453 statistic, _, _ = binned_statistic_dd(
454 [focalPlane_x[detectorInd], focalPlane_y[detectorInd]],
455 data["z"][detectorInd],
456 statistic=self.statistic,
457 bins=[1, 1],
458 )
459 patches.append(Polygon(corners, closed=True))
460 values.append(statistic.ravel()[0])
461 else:
462 # It's at amplifier level. This uses the focal
463 # plane position of the corners of the detector to
464 # generate corners for the individual amplifier
465 # segments.
466 rotation = detector.getOrientation().getNQuarter() # N * 90 degrees.
467 alpha, beta = np.cos(rotation * np.pi / 2.0), np.sin(rotation * np.pi / 2.0)
469 # Calculate the rotation matrix between X/Y and U/V
470 # coordinates.
471 scaleUX = alpha * (maxU - minU) / (maxX - minX)
472 scaleVX = beta * (maxV - minV) / (maxX - minX)
473 scaleVY = alpha * (maxV - minV) / (maxY - minY)
474 scaleUY = beta * (maxU - minU) / (maxY - minY)
476 # After the rotation, some of the corners may have
477 # negative offsets. This corresponds to corners that
478 # reference the maximum edges of the box in U/V
479 # coordinates.
480 baseU = minU if rotation % 4 in (0, 1) else maxU
481 baseV = maxV if rotation % 4 in (2, 3) else minV
483 for amplifier in detector:
484 ampName = amplifier.getName()
485 detectorInd = data["detector"] == detectorId
486 ampInd = data["amplifier"] == ampName
487 ampInd &= detectorInd
489 # Determine amplifier extent in X/Y coordinates.
490 ampMinX, ampMinY = amplifier.getBBox().getMin()
491 ampMaxX, ampMaxY = amplifier.getBBox().getMax()
493 # The corners are rotated into U/V coordinates,
494 # and the appropriate offset added.
495 ampCorners = []
496 ampCorners.append(
497 (
498 scaleUX * (ampMinX - minX) + scaleUY * (ampMinY - minY) + baseU,
499 scaleVY * (ampMinY - minY) + scaleVX * (ampMinX - minX) + baseV,
500 )
501 )
502 ampCorners.append(
503 (
504 scaleUX * (ampMaxX - minX) + scaleUY * (ampMaxY - minY) + baseU,
505 scaleVY * (ampMinY - minY) + scaleVX * (ampMinX - minX) + baseV,
506 )
507 )
508 ampCorners.append(
509 (
510 scaleUX * (ampMaxX - minX) + scaleUY * (ampMaxY - minY) + baseU,
511 scaleVY * (ampMaxY - minY) + scaleVX * (ampMaxX - minX) + baseV,
512 )
513 )
514 ampCorners.append(
515 (
516 scaleUX * (ampMinX - minX) + scaleUY * (ampMinY - minY) + baseU,
517 scaleVY * (ampMaxY - minY) + scaleVX * (ampMaxX - minX) + baseV,
518 )
519 )
520 patches.append(Polygon(ampCorners, closed=True))
521 # This does the appropriate statistic for this
522 # amplifier's data.
523 if len(data["z"][ampInd]) > 0:
524 statistic, _, _ = binned_statistic_dd(
525 [focalPlane_x[ampInd], focalPlane_y[ampInd]],
526 data["z"][ampInd],
527 statistic=self.statistic,
528 bins=[1, 1],
529 )
530 values.append(statistic.ravel()[0])
531 else:
532 values.append(np.nan)
534 # Set bounding box for this figure.
535 ax.set_xlim(plotLimit_x)
536 ax.set_ylim(plotLimit_y)
538 # Do not mask values.
539 if self.showStats:
540 statMed, statMad, statsText = self.statsAndText(values, mask=None)
541 bbox = dict(facecolor="paleturquoise", alpha=0.5, edgecolor="none")
542 ax.text(0.8, 0.91, statsText, transform=fig.transFigure, fontsize=8, bbox=bbox)
544 # Defaults to med + 4 sigma Mad to match
545 # the camera team plots
546 if self.plotMin is not None:
547 vmin = self.plotMin
548 else:
549 vmin = statMed - 4.0 * statMad
550 if self.plotMax is not None:
551 vmax = self.plotMax
552 else:
553 vmax = statMed + 4.0 * statMad
555 valuesPlot = np.clip(values, vmin, vmax)
557 patchCollection = PatchCollection(
558 patches, edgecolor="white", cmap=cmap, linewidth=0.5, linestyle=(0, (0.5, 3))
559 )
560 patchCollection.set_array(valuesPlot)
561 ax.add_collection(patchCollection)
563 divider = make_axes_locatable(ax)
564 cax = divider.append_axes("right", size="5%", pad=0.05)
565 fig.colorbar(patchCollection, cax=cax, extend="both")
566 text = cax.text(
567 0.5,
568 0.5,
569 self.zAxisLabel,
570 color="k",
571 rotation="vertical",
572 transform=cax.transAxes,
573 ha="center",
574 va="center",
575 fontsize=10,
576 )
577 text.set_path_effects([pathEffects.Stroke(linewidth=3, foreground="w"), pathEffects.Normal()])
578 cax.tick_params(labelsize=7)
580 ax.set_xlabel(self.xAxisLabel)
581 ax.set_ylabel(self.yAxisLabel)
582 ax.tick_params(axis="x", labelrotation=25)
583 ax.tick_params(labelsize=7)
585 ax.set_aspect("equal")
587 if self.addHistogram:
588 self._addHistogram(histAx, data["z"])
590 plt.draw()
592 # Add useful information to the plot
593 fig.subplots_adjust(left=0.05, right=0.95)
594 fig = plt.gcf()
595 if plotInfo:
596 fig = addPlotInfo(fig, plotInfo)
598 return fig