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