Coverage for python/lsst/analysis/tools/actions/plot/wholeSkyPlot.py: 16%
154 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-14 02:48 -0700
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__ = ("WholeSkyPlot",)
26from typing import Mapping, Optional
28import matplotlib.patheffects as pathEffects
29import matplotlib.pyplot as plt
30import numpy as np
31from lsst.pex.config import ChoiceField, Field, ListField
32from matplotlib.collections import PatchCollection
33from matplotlib.colors import CenteredNorm
34from matplotlib.figure import Figure
35from matplotlib.patches import Patch, Polygon
37from ...interfaces import KeyedData, KeyedDataSchema, PlotAction, Scalar, Vector
38from ...math import nanSigmaMad
39from ...utils import getTractCorners
40from .plotUtils import addPlotInfo, mkColormap
43class WholeSkyPlot(PlotAction):
44 """Plots the on sky distribution of a parameter.
46 Plots the values of the parameter given for the z axis
47 according to the positions given for x and y. Optimised
48 for use with RA and Dec. Also calculates some basic
49 statistics and includes those on the plot.
51 The default axes limits and figure size were chosen to plot HSC PDR2.
52 """
54 xAxisLabel = Field[str](doc="Label to use for the x axis.", default="RA (degrees)")
55 yAxisLabel = Field[str](doc="Label to use for the y axis.", default="Dec (degrees)")
56 zAxisLabel = Field[str](doc="Label to use for the z axis.", optional=False)
57 autoAxesLimits = Field[bool](doc="Find axes limits automatically.", default=True)
58 xLimits = ListField[float](doc="Plotting limits for the x axis.", default=[-5.0, 365.0])
59 yLimits = ListField[float](doc="Plotting limits for the y axis.", default=[-10.0, 60.0])
60 figureSize = ListField[float](doc="Size of the figure.", default=[9.0, 3.5])
61 colorBarRange = Field[float](
62 doc="The multiplier for the color bar range. The max/min range values are: median +/- N * sigmaMad"
63 ", where N is this config value.",
64 default=3.0,
65 )
66 colorMapType = ChoiceField[str](
67 doc="Type of color map to use for the color bar. Options: sequential, divergent, userDefined.",
68 allowed={cmType: cmType for cmType in ("sequential", "divergent")},
69 default="divergent",
70 )
71 colorMap = ListField[str](
72 doc="List of hexidecimal colors for a user-defined color map.",
73 optional=True,
74 )
76 def getInputSchema(self, **kwargs) -> KeyedDataSchema:
77 base = []
78 base.append(("z", Vector))
79 base.append(("tract", Vector))
80 return base
82 def __call__(self, data: KeyedData, **kwargs) -> Mapping[str, Figure] | Figure:
83 self._validateInput(data, **kwargs)
84 return self.makePlot(data, **kwargs)
86 def _validateInput(self, data: KeyedData, **kwargs) -> None:
87 """NOTE currently can only check that something is not a Scalar, not
88 check that the data is consistent with Vector
89 """
90 needed = self.getInputSchema(**kwargs)
91 if remainder := {key.format(**kwargs) for key, _ in needed} - {
92 key.format(**kwargs) for key in data.keys()
93 }:
94 raise ValueError(f"Task needs keys {remainder} but they were not found in input")
95 for name, typ in needed:
96 isScalar = issubclass((colType := type(data[name.format(**kwargs)])), Scalar)
97 if isScalar and typ != Scalar:
98 raise ValueError(f"Data keyed by {name} has type {colType} but action requires type {typ}")
100 def _getAxesLimits(self, xs: list, ys: list) -> tuple(list, list):
101 """Get the x and y axes limits in degrees.
103 Parameters
104 ----------
105 xs : `list`
106 X coordinates for the tracts to plot.
107 ys : `list`
108 Y coordinates for the tracts to plot.
110 Returns
111 -------
112 xlim : `list`
113 Minimun and maximum x axis values.
114 ylim : `list`
115 Minimun and maximum y axis values.
116 """
118 # Add some blank space on the edges of the plot.
119 xlim = [np.nanmin(xs) - 5, np.nanmax(xs) + 5]
120 ylim = [np.nanmin(ys) - 5, np.nanmax(ys) + 5]
122 # Limit to only show real RA/Dec values.
123 if xlim[0] < 0.0:
124 xlim[0] = 0.0
125 if xlim[1] > 360.0:
126 xlim[1] = 360.0
127 if ylim[0] < -90.0:
128 ylim[0] = -90.0
129 if ylim[1] > 90.0:
130 ylim[1] = 90.0
132 return (xlim, ylim)
134 def _getMaxOutlierVals(self, multiplier: float, tracts: list, values: list, outlierInds: list) -> str:
135 """Get the 5 largest outlier values in a string.
137 Parameters
138 ----------
139 multiplier : `float`
140 Select values whose absolute value is > multiplier * sigmaMAD.
141 tracts : `list`
142 All the tracts.
143 values : `list`
144 All the metric values.
145 outlierInds : `list`
146 Indicies of outlier values.
148 Returns
149 -------
150 text : `str`
151 A string containing the 5 tracts with the largest outlier values.
152 """
153 text = f"Tracts with |value| > {multiplier}" + r"$\sigma_{MAD}$" + ": "
154 if len(outlierInds) > 0:
155 outlierValues = np.array(values)[outlierInds]
156 outlierTracts = np.array(tracts)[outlierInds]
157 # Sort values in descending (-) absolute value order discounting
158 # NaNs.
159 maxInds = np.argsort(-np.abs(outlierValues))
160 # Show up to five values on the plot.
161 for ind in maxInds[:5]:
162 val = outlierValues[ind]
163 tract = outlierTracts[ind]
164 text += f"{tract}, {val:.3}; "
165 # Remove the final trailing comma and whitespace.
166 text = text[:-2]
167 else:
168 text += "None"
170 return text
172 def makePlot(
173 self,
174 data: KeyedData,
175 plotInfo: Optional[Mapping[str, str]] = None,
176 **kwargs,
177 ) -> Figure:
178 """Make a WholeSkyPlot of the given data.
180 Parameters
181 ----------
182 data : `KeyedData`
183 The catalog to plot the points from.
184 plotInfo : `dict`
185 A dictionary of information about the data being plotted with keys:
187 ``"run"``
188 The output run for the plots (`str`).
189 ``"skymap"``
190 The type of skymap used for the data (`str`).
191 ``"filter"``
192 The filter used for this data (`str`).
193 ``"tract"``
194 The tract that the data comes from (`str`).
196 Returns
197 -------
198 `pipeBase.Struct` containing:
199 skyPlot : `matplotlib.figure.Figure`
200 The resulting figure.
203 Examples
204 --------
205 An example of the plot produced from this code is here:
207 .. image:: /_static/analysis_tools/wholeSkyPlotExample.png
209 For a detailed example of how to make a plot from the command line
210 please see the
211 :ref:`getting started guide<analysis-tools-getting-started>`.
212 """
213 skymap = kwargs["skymap"]
214 if plotInfo is None:
215 plotInfo = {}
217 # Prevent Bands in the plot info showing a list of bands.
218 # If bands is a list, it implies that parameterizedBand=False,
219 # and that the metric is not band-specific.
220 if "bands" in plotInfo:
221 if isinstance(plotInfo["bands"], list):
222 plotInfo["bands"] = "N/A"
224 colorMap = self.colorMap
225 match self.colorMapType:
226 case "sequential":
227 if colorMap is None:
228 colorMap = mkColormap(["#F5F5F5", "#5AB4AC", "#284D48"])
229 norm = None
230 case "divergent":
231 if colorMap is None:
232 colorMap = mkColormap(["#9A6E3A", "#C6A267", "#A9A9A9", "#4F938B", "#2C665A"])
233 norm = CenteredNorm()
235 # Create patches using the corners of each tract.
236 patches = []
237 colBarVals = []
238 tracts = []
239 ras = []
240 decs = []
241 for i, tract in enumerate(data["tract"]):
242 corners = getTractCorners(skymap, tract)
243 patches.append(Polygon(corners, closed=True))
244 colBarVals.append(data["z"][i])
245 tracts.append(tract)
246 ras.append(corners[0][0])
247 decs.append(corners[0][1])
249 # Setup figure.
250 fig, ax = plt.subplots(1, 1, figsize=self.figureSize, dpi=500)
251 if self.autoAxesLimits:
252 xlim, ylim = self._getAxesLimits(ras, decs)
253 else:
254 xlim, ylim = self.xLimits, self.yLimits
255 ax.set_xlim(xlim)
256 ax.set_ylim(ylim)
257 ax.set_xlabel(self.xAxisLabel)
258 ax.set_ylabel(self.yAxisLabel)
259 ax.invert_xaxis()
261 # Add colored patches showing tract metric values.
262 patchCollection = PatchCollection(patches, cmap=colorMap, norm=norm)
263 ax.add_collection(patchCollection)
265 # Define color bar range.
266 med = np.nanmedian(colBarVals)
267 sigmaMad = nanSigmaMad(colBarVals)
268 vmin = med - self.colorBarRange * sigmaMad
269 vmax = med + self.colorBarRange * sigmaMad
271 # Note tracts with metrics outside (vmin, vmax) as outliers.
272 outlierInds = np.where((colBarVals < vmin) | (colBarVals > vmax))[0]
274 # Initialize legend handles.
275 handles = []
277 # Plot the outlier patches.
278 outlierPatches = []
279 if len(outlierInds) > 0:
280 for ind in outlierInds:
281 outlierPatches.append(patches[ind])
282 outlierPatchCollection = PatchCollection(
283 outlierPatches,
284 cmap=colorMap,
285 norm=norm,
286 facecolors="none",
287 edgecolors="k",
288 linewidths=0.5,
289 zorder=100,
290 )
291 ax.add_collection(outlierPatchCollection)
292 # Add legend information.
293 outlierPatch = Patch(
294 facecolor="none",
295 edgecolor="k",
296 linewidth=0.5,
297 label="Outlier",
298 )
299 handles.append(outlierPatch)
301 # Plot tracts with NaN metric values.
302 nanInds = np.where(~np.isfinite(colBarVals))[0]
303 nanPatches = []
304 if len(nanInds) > 0:
305 for ind in nanInds:
306 nanPatches.append(patches[ind])
307 nanPatchCollection = PatchCollection(
308 nanPatches,
309 cmap=None,
310 norm=norm,
311 facecolors="white",
312 edgecolors="grey",
313 linestyles="dotted",
314 linewidths=0.5,
315 zorder=10,
316 )
317 ax.add_collection(nanPatchCollection)
318 # Add legend information.
319 nanPatch = Patch(
320 facecolor="white",
321 edgecolor="grey",
322 linestyle="dotted",
323 linewidth=0.5,
324 label="NaN",
325 )
326 handles.append(nanPatch)
328 if len(handles) > 0:
329 fig.legend(handles=handles)
331 # Add text boxes to show the number of tracts, number of NaNs,
332 # median, sigma MAD, and the five largest outlier values.
333 outlierText = self._getMaxOutlierVals(self.colorBarRange, tracts, colBarVals, outlierInds)
334 # Make vertical text spacing readable for different figure sizes.
335 multiplier = 3.5 / self.figureSize[1]
336 verticalSpacing = 0.028 * multiplier
337 fig.text(
338 0.01,
339 0.01 + 3 * verticalSpacing,
340 f"Num tracts: {len(tracts)}",
341 transform=fig.transFigure,
342 fontsize=8,
343 alpha=0.7,
344 )
345 fig.text(
346 0.01,
347 0.01 + 2 * verticalSpacing,
348 f"Num nans: {len(nanInds)}",
349 transform=fig.transFigure,
350 fontsize=8,
351 alpha=0.7,
352 )
353 fig.text(
354 0.01,
355 0.01 + verticalSpacing,
356 f"Median: {med:.3f}; " + r"$\sigma_{MAD}$" + f": {sigmaMad:.3f}",
357 transform=fig.transFigure,
358 fontsize=8,
359 alpha=0.7,
360 )
361 fig.text(0.01, 0.01, outlierText, transform=fig.transFigure, fontsize=8, alpha=0.7)
363 # Truncate the color range to (vmin, vmax).
364 colorBarVals = np.clip(np.array(colBarVals), vmin, vmax)
365 patchCollection.set_array(colorBarVals)
366 # Make the color bar with a metric label.
367 cbar = plt.colorbar(
368 patchCollection,
369 ax=ax,
370 shrink=0.7,
371 extend="both",
372 location="top",
373 orientation="horizontal",
374 )
375 cbarText = self.zAxisLabel.format_map(kwargs)
376 if "zUnit" in data and data["zUnit"] != "":
377 cbarText += f" ({data['zUnit']})"
378 text = cbar.ax.text(
379 0.5,
380 0.5,
381 cbarText,
382 transform=cbar.ax.transAxes,
383 ha="center",
384 va="center",
385 fontsize=10,
386 zorder=100,
387 )
388 text.set_path_effects([pathEffects.Stroke(linewidth=3, foreground="w"), pathEffects.Normal()])
390 # Finalize plot appearance.
391 ax.grid()
392 ax.set_axisbelow(True)
393 ax.set_aspect("equal")
394 fig = plt.gcf()
395 fig = addPlotInfo(fig, plotInfo)
396 plt.subplots_adjust(left=0.08, right=0.97, top=0.8, bottom=0.17, wspace=0.35)
398 return fig