From 65dccfc4501ea8400c5afb23a2da16a87c3d7be4 Mon Sep 17 00:00:00 2001 From: Christopher Waters Date: Fri, 26 May 2023 17:02:04 -0700 Subject: [PATCH 1/5] Add class to plot per-detector or per-amplifier values for the full focal plane. --- .../tools/actions/plot/focalPlanePlot.py | 240 +++++++++++++++++- 1 file changed, 237 insertions(+), 3 deletions(-) diff --git a/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py b/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py index 4d7ad3e59..4a2b68aa1 100644 --- a/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py +++ b/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py @@ -21,7 +21,7 @@ from __future__ import annotations -__all__ = ("FocalPlanePlot",) +__all__ = ("FocalPlanePlot", "FocalPlaneGeometryPlot") from typing import Mapping, Optional @@ -29,9 +29,11 @@ import matplotlib.pyplot as plt import numpy as np from lsst.afw.cameraGeom import FOCAL_PLANE, PIXELS, Camera -from lsst.pex.config import Field +from lsst.pex.config import ChoiceField, Field +from matplotlib.collections import PatchCollection from matplotlib.figure import Figure -from scipy.stats import binned_statistic_2d +from matplotlib.patches import Polygon +from scipy.stats import binned_statistic_2d, binned_statistic_dd from ...interfaces import KeyedData, KeyedDataSchema, PlotAction, Scalar, Vector from ...statistics import nansigmaMad @@ -230,3 +232,235 @@ def makePlot( fig = addPlotInfo(fig, plotInfo) return fig + + +class FocalPlaneGeometryPlot(FocalPlanePlot): + """Plots the focal plane distribution of a parameter in camera + geometry units. + + Given the detector positions in x and y, the focal plane positions + are calculated using the camera model. A 2d binned statistic + (default is mean) is then calculated and plotted for the parameter + z as a function of the camera geometry segment the input points + fall upon. + """ + + xAxisLabel = Field[str](doc="Label to use for the x axis.", optional=False) + yAxisLabel = Field[str](doc="Label to use for the y axis.", optional=False) + zAxisLabel = Field[str](doc="Label to use for the z axis.", optional=False) + + level = ChoiceField[str]( + doc="Which geometry level should values be plotted?", + default="amplifier", + allowed={ + "amplifier": "Plot values per readout amplifier.", + "detector": "Plot values per detector.", + }, + ) + statistic = Field[str]( + doc="Operation to perform in binned_statistic_2d", + default="mean", + ) + + def makePlot( + self, + data: KeyedData, + camera: Camera, + plotInfo: Optional[Mapping[str, str]] = None, + **kwargs, + ) -> Figure: + """Prep the catalogue and then make a focalPlanePlot of the given + column. + + Uses the axisLabels config options `x` and `y` to make an image, where + the color corresponds to the 2d binned statistic (the mean is the + default) applied to the `z` column. A summary panel is shown in the + upper right corner of the resultant plot. The code uses the + selectorActions to decide which points to plot and the + statisticSelector actions to determine which points to use for the + printed statistics. + + Parameters + ---------- + data : `pandas.core.frame.DataFrame` + The catalog to plot the points from. + camera : `lsst.afw.cameraGeom.Camera` + The camera used to map from pixel to focal plane positions. + plotInfo : `dict` + A dictionary of information about the data being plotted with keys: + ``"run"`` + The output run for the plots (`str`). + ``"skymap"`` + The type of skymap used for the data (`str`). + ``"filter"`` + The filter used for this data (`str`). + ``"tract"`` + The tract that the data comes from (`str`). + ``"bands"`` + The band(s) that the data comes from (`list` of `str`). + + Returns + ------- + fig : `matplotlib.figure.Figure` + The resulting figure. + """ + if plotInfo is None: + plotInfo = {} + + if len(data["x"]) == 0: + noDataFig = Figure() + noDataFig.text(0.3, 0.5, "No data to plot after selectors applied") + noDataFig = addPlotInfo(noDataFig, plotInfo) + return noDataFig + + fig = plt.figure(dpi=300) + ax = fig.add_subplot(111) + + detectorIds = np.unique(data["detector"]) + focalPlane_x = np.zeros(len(data["x"])) + focalPlane_y = np.zeros(len(data["y"])) + + patches = [] + values = [] + + # Plot bounding box that will be used to set the axes below. + plotLimit_x = [0.0, 0.0] + plotLimit_y = [0.0, 0.0] + + for detectorId in detectorIds: + detector = camera[detectorId] + # We can go stright to fp coordinates. + corners = [(c.getX(), c.getY()) for c in detector.getCorners(FOCAL_PLANE)] + corners = np.array(corners) + + # See if the plot bounding box needs to be extended: + for corner in corners: + if corner[0] < plotLimit_x[0]: + plotLimit_x[0] = corner[0] + if corner[0] > plotLimit_x[1]: + plotLimit_x[1] = corner[0] + if corner[1] < plotLimit_y[0]: + plotLimit_y[0] = corner[1] + if corner[1] > plotLimit_y[1]: + plotLimit_y[1] = corner[1] + + # U/V coordinates represent focal plane locations. + minU, minV = corners.min(axis=0) + maxU, maxV = corners.max(axis=0) + + # X/Y coordinates represent detector internal coordinates. + # Detector extent in detector coordinates + minX, minY = detector.getBBox().getMin() + maxX, maxY = detector.getBBox().getMax() + + if self.level.lower() == 'detector': + detectorInd = data["detector"] == detectorId + + # This does the appropriate statistic for this + # detector's data. + statistic, _, _ = binned_statistic_dd([focalPlane_x[detectorInd], + focalPlane_y[detectorInd]], + data["z"][detectorInd], + statistic=self.statistic, + bins=[1, 1]) + patches.append(Polygon(corners, True)) + values.append(statistic.ravel()[0]) + else: + # It's at amplifier level. This uses the focal + # plane position of the corners of the detector to + # generate corners for the individual amplifier + # segments. + rotation = detector.getOrientation().getNQuarter() # N * 90 degrees. + alpha, beta = np.cos(rotation * np.pi / 2.0), np.sin(rotation * np.pi / 2.0) + + # Calculate the rotation matrix between X/Y and U/V + # coordinates. + scaleUX = alpha * (maxU - minU)/(maxX - minX) + scaleVX = beta * (maxV - minV)/(maxX - minX) + scaleVY = alpha * (maxV - minV)/(maxY - minY) + scaleUY = beta * (maxU - minU)/(maxY - minY) + + # After the rotation, some of the corners may have + # negative offsets. This corresponds to corners that + # reference the maximum edges of the box in U/V + # coordinates. + baseU = minU if rotation % 4 in (0, 1) else maxU + baseV = maxV if rotation % 4 in (2, 3) else minV + + for amplifier in detector: + ampName = amplifier.getName() + detectorInd = data["detector"] == detectorId + ampInd = data["x"] == ampName + detectorInd &= ampInd + + # Determine amplifier extent in X/Y coordinates. + ampMinX, ampMinY = amplifier.getBBox().getMin() + ampMaxX, ampMaxY = amplifier.getBBox().getMax() + + # The corners are rotated into U/V coordinates, + # and the appropriate offset added. + ampCorners = [] + ampCorners.append((scaleUX * (ampMinX - minX) + scaleUY * (ampMinY - minY) + baseU, + scaleVY * (ampMinY - minY) + scaleVX * (ampMinX - minX) + baseV)) + ampCorners.append((scaleUX * (ampMaxX - minX) + scaleUY * (ampMaxY - minY) + baseU, + scaleVY * (ampMinY - minY) + scaleVX * (ampMinX - minX) + baseV)) + ampCorners.append((scaleUX * (ampMaxX - minX) + scaleUY * (ampMaxY - minY) + baseU, + scaleVY * (ampMaxY - minY) + scaleVX * (ampMaxX - minX) + baseV)) + ampCorners.append((scaleUX * (ampMinX - minX) + scaleUY * (ampMinY - minY) + baseU, + scaleVY * (ampMaxY - minY) + scaleVX * (ampMaxX - minX) + baseV)) + patches.append(Polygon(ampCorners, True)) + # This does the appropriate statistic for this + # amplifier's data. + if len(data["z"][detectorInd]) > 0: + statistic, _, _ = binned_statistic_dd([focalPlane_x[detectorInd], + focalPlane_y[detectorInd]], + data["z"][detectorInd], + statistic=self.statistic, + bins=[1, 1]) + values.append(statistic.ravel()[0]) + else: + values.append(np.nan) + + # Set bounding box for this figure. + ax.set_xlim(plotLimit_x) + ax.set_ylim(plotLimit_y) + + # Do not mask values. + statMed, statMad, statsText = self.statsAndText(values, mask=None) + bbox = dict(facecolor="paleturquoise", alpha=0.5, edgecolor="none") + ax.text(0.8, 0.91, statsText, transform=fig.transFigure, fontsize=8, bbox=bbox) + + patchCollection = PatchCollection(patches, alpha=0.4, edgecolor='black') + patchCollection.set_array(values) + ax.add_collection(patchCollection) + + cax = fig.add_axes([0.87 + 0.04, 0.11, 0.04, 0.77]) + fig.colorbar(patchCollection, cax=cax, extend="both") + text = cax.text( + 0.5, + 0.5, + self.zAxisLabel, + color="k", + rotation="vertical", + transform=cax.transAxes, + ha="center", + va="center", + fontsize=10, + ) + text.set_path_effects([pathEffects.Stroke(linewidth=3, foreground="w"), pathEffects.Normal()]) + cax.tick_params(labelsize=7) + + ax.set_xlabel(self.xAxisLabel) + ax.set_ylabel(self.yAxisLabel) + ax.tick_params(axis="x", labelrotation=25) + ax.tick_params(labelsize=7) + + ax.set_aspect("equal") + plt.draw() + + # Add useful information to the plot + plt.subplots_adjust(wspace=0.0, hspace=0.0, right=0.85) + fig = plt.gcf() + fig = addPlotInfo(fig, plotInfo) + + return fig From 93481749deaaa6df2b226423a98be8deb8d693e7 Mon Sep 17 00:00:00 2001 From: Christopher Waters Date: Mon, 10 Jul 2023 15:01:52 -0700 Subject: [PATCH 2/5] Fix formatting. --- .../tools/actions/plot/focalPlanePlot.py | 66 ++++++++++++------- 1 file changed, 42 insertions(+), 24 deletions(-) diff --git a/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py b/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py index 4a2b68aa1..d430968fe 100644 --- a/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py +++ b/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py @@ -353,16 +353,17 @@ def makePlot( minX, minY = detector.getBBox().getMin() maxX, maxY = detector.getBBox().getMax() - if self.level.lower() == 'detector': + if self.level.lower() == "detector": detectorInd = data["detector"] == detectorId # This does the appropriate statistic for this # detector's data. - statistic, _, _ = binned_statistic_dd([focalPlane_x[detectorInd], - focalPlane_y[detectorInd]], - data["z"][detectorInd], - statistic=self.statistic, - bins=[1, 1]) + statistic, _, _ = binned_statistic_dd( + [focalPlane_x[detectorInd], focalPlane_y[detectorInd]], + data["z"][detectorInd], + statistic=self.statistic, + bins=[1, 1], + ) patches.append(Polygon(corners, True)) values.append(statistic.ravel()[0]) else: @@ -375,10 +376,10 @@ def makePlot( # Calculate the rotation matrix between X/Y and U/V # coordinates. - scaleUX = alpha * (maxU - minU)/(maxX - minX) - scaleVX = beta * (maxV - minV)/(maxX - minX) - scaleVY = alpha * (maxV - minV)/(maxY - minY) - scaleUY = beta * (maxU - minU)/(maxY - minY) + scaleUX = alpha * (maxU - minU) / (maxX - minX) + scaleVX = beta * (maxV - minV) / (maxX - minX) + scaleVY = alpha * (maxV - minV) / (maxY - minY) + scaleUY = beta * (maxU - minU) / (maxY - minY) # After the rotation, some of the corners may have # negative offsets. This corresponds to corners that @@ -400,23 +401,40 @@ def makePlot( # The corners are rotated into U/V coordinates, # and the appropriate offset added. ampCorners = [] - ampCorners.append((scaleUX * (ampMinX - minX) + scaleUY * (ampMinY - minY) + baseU, - scaleVY * (ampMinY - minY) + scaleVX * (ampMinX - minX) + baseV)) - ampCorners.append((scaleUX * (ampMaxX - minX) + scaleUY * (ampMaxY - minY) + baseU, - scaleVY * (ampMinY - minY) + scaleVX * (ampMinX - minX) + baseV)) - ampCorners.append((scaleUX * (ampMaxX - minX) + scaleUY * (ampMaxY - minY) + baseU, - scaleVY * (ampMaxY - minY) + scaleVX * (ampMaxX - minX) + baseV)) - ampCorners.append((scaleUX * (ampMinX - minX) + scaleUY * (ampMinY - minY) + baseU, - scaleVY * (ampMaxY - minY) + scaleVX * (ampMaxX - minX) + baseV)) + ampCorners.append( + ( + scaleUX * (ampMinX - minX) + scaleUY * (ampMinY - minY) + baseU, + scaleVY * (ampMinY - minY) + scaleVX * (ampMinX - minX) + baseV, + ) + ) + ampCorners.append( + ( + scaleUX * (ampMaxX - minX) + scaleUY * (ampMaxY - minY) + baseU, + scaleVY * (ampMinY - minY) + scaleVX * (ampMinX - minX) + baseV, + ) + ) + ampCorners.append( + ( + scaleUX * (ampMaxX - minX) + scaleUY * (ampMaxY - minY) + baseU, + scaleVY * (ampMaxY - minY) + scaleVX * (ampMaxX - minX) + baseV, + ) + ) + ampCorners.append( + ( + scaleUX * (ampMinX - minX) + scaleUY * (ampMinY - minY) + baseU, + scaleVY * (ampMaxY - minY) + scaleVX * (ampMaxX - minX) + baseV, + ) + ) patches.append(Polygon(ampCorners, True)) # This does the appropriate statistic for this # amplifier's data. if len(data["z"][detectorInd]) > 0: - statistic, _, _ = binned_statistic_dd([focalPlane_x[detectorInd], - focalPlane_y[detectorInd]], - data["z"][detectorInd], - statistic=self.statistic, - bins=[1, 1]) + statistic, _, _ = binned_statistic_dd( + [focalPlane_x[detectorInd], focalPlane_y[detectorInd]], + data["z"][detectorInd], + statistic=self.statistic, + bins=[1, 1], + ) values.append(statistic.ravel()[0]) else: values.append(np.nan) @@ -430,7 +448,7 @@ def makePlot( bbox = dict(facecolor="paleturquoise", alpha=0.5, edgecolor="none") ax.text(0.8, 0.91, statsText, transform=fig.transFigure, fontsize=8, bbox=bbox) - patchCollection = PatchCollection(patches, alpha=0.4, edgecolor='black') + patchCollection = PatchCollection(patches, alpha=0.4, edgecolor="black") patchCollection.set_array(values) ax.add_collection(patchCollection) From f71cff7230e17c0c783715e67af1055250e83ff5 Mon Sep 17 00:00:00 2001 From: Christopher Waters Date: Fri, 21 Jul 2023 14:26:31 -0700 Subject: [PATCH 3/5] Update docstrings. --- .../tools/actions/plot/focalPlanePlot.py | 29 ++++++++++++------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py b/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py index d430968fe..9a5a3d2c0 100644 --- a/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py +++ b/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py @@ -235,19 +235,18 @@ def makePlot( class FocalPlaneGeometryPlot(FocalPlanePlot): - """Plots the focal plane distribution of a parameter in camera - geometry units. + """Plots the focal plane distribution of a parameter in afw camera + geometry units: amplifiers and detectors. Given the detector positions in x and y, the focal plane positions are calculated using the camera model. A 2d binned statistic (default is mean) is then calculated and plotted for the parameter z as a function of the camera geometry segment the input points fall upon. - """ - xAxisLabel = Field[str](doc="Label to use for the x axis.", optional=False) - yAxisLabel = Field[str](doc="Label to use for the y axis.", optional=False) - zAxisLabel = Field[str](doc="Label to use for the z axis.", optional=False) + The ``xAxisLabel``, ``yAxisLabel``, ``zAxisLabel``, and + ``statistic`` variables are inherited from the parent class. + """ level = ChoiceField[str]( doc="Which geometry level should values be plotted?", @@ -257,10 +256,6 @@ class FocalPlaneGeometryPlot(FocalPlanePlot): "detector": "Plot values per detector.", }, ) - statistic = Field[str]( - doc="Operation to perform in binned_statistic_2d", - default="mean", - ) def makePlot( self, @@ -283,7 +278,19 @@ def makePlot( Parameters ---------- data : `pandas.core.frame.DataFrame` - The catalog to plot the points from. + The catalog to plot the points from. This is expected to + have the following columns/keys: + ``"detector"`` + The integer detector id for the points. + ``"amplifier"`` + The string amplifier name for the points. + ``"z"`` + The numerical value that will be combined via + ``statistic`` to the binned value. + ``"x"`` + Focal plane x position, optional. + ``"y"`` + Focal plane y position, optional. camera : `lsst.afw.cameraGeom.Camera` The camera used to map from pixel to focal plane positions. plotInfo : `dict` From d5a8ca09cba4deae8624abd93da8fb20cf7eff10 Mon Sep 17 00:00:00 2001 From: Christopher Waters Date: Fri, 21 Jul 2023 14:28:09 -0700 Subject: [PATCH 4/5] Simplify and correct some variable/column names. * Switch to an explicit 'amplifier' column, so as not to overload 'x'. * Switch to using 'z' column for all expected lengths, as it should always exist. * Clarify index name, as it applies to amplifiers, not detectors. --- .../tools/actions/plot/focalPlanePlot.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py b/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py index 9a5a3d2c0..03a1af74f 100644 --- a/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py +++ b/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py @@ -314,7 +314,7 @@ def makePlot( if plotInfo is None: plotInfo = {} - if len(data["x"]) == 0: + if len(data["z"]) == 0: noDataFig = Figure() noDataFig.text(0.3, 0.5, "No data to plot after selectors applied") noDataFig = addPlotInfo(noDataFig, plotInfo) @@ -324,8 +324,8 @@ def makePlot( ax = fig.add_subplot(111) detectorIds = np.unique(data["detector"]) - focalPlane_x = np.zeros(len(data["x"])) - focalPlane_y = np.zeros(len(data["y"])) + focalPlane_x = np.zeros(len(data["z"])) + focalPlane_y = np.zeros(len(data["z"])) patches = [] values = [] @@ -398,8 +398,8 @@ def makePlot( for amplifier in detector: ampName = amplifier.getName() detectorInd = data["detector"] == detectorId - ampInd = data["x"] == ampName - detectorInd &= ampInd + ampInd = data["amplifier"] == ampName + ampInd &= detectorInd # Determine amplifier extent in X/Y coordinates. ampMinX, ampMinY = amplifier.getBBox().getMin() @@ -435,10 +435,10 @@ def makePlot( patches.append(Polygon(ampCorners, True)) # This does the appropriate statistic for this # amplifier's data. - if len(data["z"][detectorInd]) > 0: + if len(data["z"][ampInd]) > 0: statistic, _, _ = binned_statistic_dd( - [focalPlane_x[detectorInd], focalPlane_y[detectorInd]], - data["z"][detectorInd], + [focalPlane_x[ampInd], focalPlane_y[ampInd]], + data["z"][ampInd], statistic=self.statistic, bins=[1, 1], ) From 04d6015e9835bc27c8e504e9707bf50f7c4edd69 Mon Sep 17 00:00:00 2001 From: Christopher Waters Date: Fri, 21 Jul 2023 14:30:47 -0700 Subject: [PATCH 5/5] Remove unneeded iteration. --- .../tools/actions/plot/focalPlanePlot.py | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py b/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py index 03a1af74f..3b3e0b712 100644 --- a/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py +++ b/python/lsst/analysis/tools/actions/plot/focalPlanePlot.py @@ -336,25 +336,25 @@ def makePlot( for detectorId in detectorIds: detector = camera[detectorId] + # We can go stright to fp coordinates. corners = [(c.getX(), c.getY()) for c in detector.getCorners(FOCAL_PLANE)] corners = np.array(corners) - # See if the plot bounding box needs to be extended: - for corner in corners: - if corner[0] < plotLimit_x[0]: - plotLimit_x[0] = corner[0] - if corner[0] > plotLimit_x[1]: - plotLimit_x[1] = corner[0] - if corner[1] < plotLimit_y[0]: - plotLimit_y[0] = corner[1] - if corner[1] > plotLimit_y[1]: - plotLimit_y[1] = corner[1] - # U/V coordinates represent focal plane locations. minU, minV = corners.min(axis=0) maxU, maxV = corners.max(axis=0) + # See if the plot bounding box needs to be extended: + if minU < plotLimit_x[0]: + plotLimit_x[0] = minU + if minV < plotLimit_y[0]: + plotLimit_y[0] = minV + if maxU > plotLimit_x[1]: + plotLimit_x[1] = maxU + if maxV > plotLimit_y[1]: + plotLimit_y[1] = maxV + # X/Y coordinates represent detector internal coordinates. # Detector extent in detector coordinates minX, minY = detector.getBBox().getMin()