Add histograms to wholeSkyPlot, bonus speed ups, metric thresholds and error catching

sr525 · sr525 · commit 8694c7b2833c · 2025-11-20T11:38:29.000-08:00
diff --git a/python/lsst/analysis/tools/actions/plot/wholeSkyPlot.py b/python/lsst/analysis/tools/actions/plot/wholeSkyPlot.py
@@ -23,12 +23,23 @@
 
 __all__ = ("WholeSkyPlot",)
 
+import importlib.resources as importResources
 from typing import Mapping, Optional
 
+import lsst.analysis.tools
 import matplotlib.patheffects as pathEffects
-import matplotlib.pyplot as plt
 import numpy as np
+import yaml
 from lsst.pex.config import ChoiceField, Field, ListField
+from lsst.utils.plotting import (
+    accent_color,
+    divergent_cmap,
+    make_figure,
+    set_rubin_plotstyle,
+    stars_cmap,
+    stars_color,
+)
+from matplotlib import gridspec
 from matplotlib.collections import PatchCollection
 from matplotlib.colors import CenteredNorm
 from matplotlib.figure import Figure
@@ -37,7 +48,7 @@
 from ...interfaces import KeyedData, KeyedDataSchema, PlotAction, Scalar, Vector
 from ...math import nanSigmaMad
 from ...utils import getTractCorners
-from .plotUtils import addPlotInfo, mkColormap
+from .plotUtils import addPlotInfo
 
 
 class WholeSkyPlot(PlotAction):
@@ -58,8 +69,6 @@ class WholeSkyPlot(PlotAction):
     xLimits = ListField[float](doc="Plotting limits for the x axis.", default=[-5.0, 365.0])
     yLimits = ListField[float](doc="Plotting limits for the y axis.", default=[-10.0, 60.0])
     autoAxesLimits = Field[bool](doc="Find axes limits automatically.", default=True)
-    dpi = Field[int](doc="DPI size of the figure.", default=500)
-    figureSize = ListField[float](doc="Size of the figure.", default=[9.0, 3.5])
     colorBarMin = Field[float](doc="The minimum value of the color bar.", optional=True)
     colorBarMax = Field[float](doc="The minimum value of the color bar.", optional=True)
     colorBarRange = Field[float](
@@ -84,6 +93,11 @@ class WholeSkyPlot(PlotAction):
     showNaNs = Field[bool](doc="Show the NaNs on the plot.", default=True)
     labelTracts = Field[bool](doc="Label the tracts.", default=False)
 
+    addThreshold = Field[bool](
+        doc="Read in the predefined threshold and add it to the histogram.",
+        default=True,
+    )
+
     def getInputSchema(self, **kwargs) -> KeyedDataSchema:
         base = []
         base.append(("z", Vector))
@@ -161,15 +175,18 @@ def _getMaxOutlierVals(self, multiplier: float, tracts: list, values: list, outl
         text : `str`
             A string containing the 5 tracts with the largest outlier values.
         """
-        text = f"Tracts with |value| > {multiplier}" + r"$\sigma_{MAD}$" + ": "
+        if self.addThreshold:
+            text = "Tracts with value outside thresholds: "
+        else:
+            text = f"Tracts with |value| > {multiplier}" + r"$\sigma_{MAD}$" + ": "
         if len(outlierInds) > 0:
             outlierValues = np.array(values)[outlierInds]
             outlierTracts = np.array(tracts)[outlierInds]
             # Sort values in descending (-) absolute value order discounting
             # NaNs.
             maxInds = np.argsort(-np.abs(outlierValues))
             # Show up to five values on the plot.
-            for ind in maxInds[:5]:
+            for ind in maxInds[:10]:
                 val = outlierValues[ind]
                 tract = outlierTracts[ind]
                 text += f"{tract}, {val:.3}; "
@@ -225,6 +242,10 @@ def makePlot(
         if plotInfo is None:
             plotInfo = {}
 
+        if self.addThreshold:
+            metricThresholdFile = importResources.read_text(lsst.analysis.tools, "metricInformation.yaml")
+            metricDefs = yaml.safe_load(metricThresholdFile)
+
         # Prevent Bands in the plot info showing a list of bands.
         # If bands is a list, it implies that parameterizedBand=False,
         # and that the metric is not band-specific.
@@ -236,11 +257,13 @@ def makePlot(
         match self.colorMapType:
             case "sequential":
                 if colorMap is None:
-                    colorMap = mkColormap(["#F5F5F5", "#5AB4AC", "#284D48"])
+                    colorMap = stars_cmap()
+                    outlierColor = "red"
                 norm = None
             case "divergent":
                 if colorMap is None:
-                    colorMap = mkColormap(["#9A6E3A", "#C6A267", "#A9A9A9", "#4F938B", "#2C665A"])
+                    colorMap = divergent_cmap()
+                    outlierColor = "fuchsia"
                 norm = CenteredNorm()
 
         # Create patches using the corners of each tract.
@@ -262,23 +285,19 @@ def makePlot(
             mid_decs.append((corners[0][1] + corners[2][1]) / 2)
 
         # Setup figure.
-        fig, ax = plt.subplots(1, 1, figsize=self.figureSize, dpi=self.dpi)
-        if self.autoAxesLimits:
-            xlim, ylim = self._getAxesLimits(ras, decs)
-        else:
-            xlim, ylim = self.xLimits, self.yLimits
-        ax.set_xlim(xlim)
-        ax.set_ylim(ylim)
-        ax.set_xlabel(self.xAxisLabel)
-        ax.set_ylabel(self.yAxisLabel)
-        ax.invert_xaxis()
-
+        fig = make_figure(dpi=300, figsize=(12, 3.5))
+        set_rubin_plotstyle()
+        gs = gridspec.GridSpec(1, 4)
+        ax = fig.add_subplot(gs[:3])
         # Add colored patches showing tract metric values.
         patchCollection = PatchCollection(patches, cmap=colorMap, norm=norm)
         ax.add_collection(patchCollection)
 
         # Define color bar range.
-        med = np.nanmedian(colBarVals)
+        if np.sum(np.isfinite(colBarVals)) > 0:
+            med = np.nanmedian(colBarVals)
+        else:
+            med = np.nan
         sigmaMad = nanSigmaMad(colBarVals)
         if self.colorBarMin is not None:
             vmin = np.float64(self.colorBarMin)
@@ -289,8 +308,15 @@ def makePlot(
         else:
             vmax = med + self.colorBarRange * sigmaMad
 
-        # Note tracts with metrics outside (vmin, vmax) as outliers.
-        outlierInds = np.where((colBarVals < vmin) | (colBarVals > vmax))[0]
+        dataName = self.zAxisLabel.format_map(kwargs)
+        colBarVals = np.array(colBarVals)
+        if self.addThreshold and dataName in metricDefs:
+            lowThreshold = metricDefs[dataName]["lowThreshold"]
+            highThreshold = metricDefs[dataName]["highThreshold"]
+            outlierInds = np.where((colBarVals < lowThreshold) | (colBarVals > highThreshold))[0]
+        else:
+            # Note tracts with metrics outside (vmin, vmax) as outliers.
+            outlierInds = np.where((colBarVals < vmin) | (colBarVals > vmax))[0]
 
         # Initialize legend handles.
         handles = []
@@ -306,15 +332,15 @@ def makePlot(
                     cmap=colorMap,
                     norm=norm,
                     facecolors="none",
-                    edgecolors="k",
+                    edgecolors=outlierColor,
                     linewidths=0.5,
                     zorder=100,
                 )
                 ax.add_collection(outlierPatchCollection)
                 # Add legend information.
                 outlierPatch = Patch(
                     facecolor="none",
-                    edgecolor="k",
+                    edgecolor=outlierColor,
                     linewidth=0.5,
                     label="Outlier",
                 )
@@ -365,13 +391,50 @@ def makePlot(
                     zorder=100,
                 )
 
+        ax.set_aspect("equals")
+        axPos = ax.get_position()
+        ax1 = fig.add_axes([0.73, 0.25, 0.20, 0.47])
+
+        if np.sum(np.isfinite(data["z"])) > 0:
+            ax1.hist(data["z"], bins=len(data["z"] / 10), color=stars_color(), histtype="step")
+        else:
+            ax1.text(0.5, 0.5, "Data all NaN/Inf")
+        ax1.set_xlabel("Metric Values")
+        ax1.set_ylabel("Number")
+        ax1.yaxis.set_label_position("right")
+        ax1.yaxis.tick_right()
+
+        if self.addThreshold and dataName in metricDefs:
+            ax1.axvline(lowThreshold, color=accent_color())
+            ax1.axvline(highThreshold, color=accent_color())
+
+            widthThreshold = highThreshold - lowThreshold
+            upperLim = highThreshold + 0.5 * widthThreshold
+            lowerLim = lowThreshold - 0.5 * widthThreshold
+            ax1.set_xlim(lowerLim, upperLim)
+            numOutside = np.sum(((data["z"] > upperLim) | (data["z"] < lowerLim)))
+            ax1.set_title("Outside plot limits: " + str(numOutside))
+
+        else:
+            if vmin != vmax and np.isfinite(vmin) and np.isfinite(vmax):
+                ax1.set_xlim(vmin, vmax)
+
+        if self.autoAxesLimits:
+            xlim, ylim = self._getAxesLimits(ras, decs)
+        else:
+            xlim, ylim = self.xLimits, self.yLimits
+        ax.set_xlim(xlim)
+        ax.set_ylim(ylim)
+        ax.set_xlabel(self.xAxisLabel)
+        ax.set_ylabel(self.yAxisLabel)
+        ax.invert_xaxis()
+
         if self.showOutliers:
             # Add text boxes to show the number of tracts, number of NaNs,
             # median, sigma MAD, and the five largest outlier values.
             outlierText = self._getMaxOutlierVals(self.colorBarRange, tracts, colBarVals, outlierInds)
-
         # Make vertical text spacing readable for different figure sizes.
-        multiplier = 3.5 / self.figureSize[1]
+        multiplier = 3.5 / fig.get_size_inches()[1]
         verticalSpacing = 0.028 * multiplier
         fig.text(
             0.01,
@@ -402,25 +465,27 @@ def makePlot(
             fig.text(0.01, 0.01, outlierText, transform=fig.transFigure, fontsize=8, alpha=0.7)
 
         # Truncate the color range to (vmin, vmax).
-        colorBarVals = np.clip(np.array(colBarVals), vmin, vmax)
-        patchCollection.set_array(colorBarVals)
+        if vmin != vmax and np.isfinite(vmin) and np.isfinite(vmax):
+            colBarVals = np.clip(np.array(colBarVals), vmin, vmax)
+        patchCollection.set_array(colBarVals)
         # Make the color bar with a metric label.
-        cbar = plt.colorbar(
+        axPos = ax.get_position()
+        cax = fig.add_axes([0.084, axPos.y1 + 0.02, 0.62, 0.07])
+        fig.colorbar(
             patchCollection,
-            ax=ax,
+            cax=cax,
             shrink=0.7,
             extend="both",
             location="top",
             orientation="horizontal",
         )
-        cbarText = self.zAxisLabel.format_map(kwargs)
-        if "zUnit" in data and data["zUnit"] != "":
-            cbarText += f" ({data['zUnit']})"
-        text = cbar.ax.text(
+        cbarText = "Metric Values"
+
+        text = cax.text(
             0.5,
             0.5,
             cbarText,
-            transform=cbar.ax.transAxes,
+            transform=cax.transAxes,
             ha="center",
             va="center",
             fontsize=10,
@@ -431,9 +496,11 @@ def makePlot(
         # Finalize plot appearance.
         ax.grid()
         ax.set_axisbelow(True)
-        ax.set_aspect("equal")
-        fig = plt.gcf()
         fig = addPlotInfo(fig, plotInfo)
-        plt.subplots_adjust(left=0.08, right=0.97, top=0.8, bottom=0.17, wspace=0.35)
+        fig.subplots_adjust(left=0.08, right=0.92, top=0.8, bottom=0.17, wspace=0.05)
+        titleText = self.zAxisLabel.format_map(kwargs)
+        if "zUnit" in data and data["zUnit"] != "":
+            titleText += f" ({data['zUnit']})"
+        fig.suptitle("Metric: " + titleText, fontsize=20)
 
         return fig
