Coverage for python/lsst/analysis/tools/actions/plot/wholeSkyPlot.py: 16%
177 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-22 01:25 -0800
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.", default="")
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 autoAxesLimits = Field[bool](doc="Find axes limits automatically.", default=True)
61 dpi = Field[int](doc="DPI size of the figure.", default=500)
62 figureSize = ListField[float](doc="Size of the figure.", default=[9.0, 3.5])
63 colorBarMin = Field[float](doc="The minimum value of the color bar.", optional=True)
64 colorBarMax = Field[float](doc="The minimum value of the color bar.", optional=True)
65 colorBarRange = Field[float](
66 doc="The multiplier for the color bar range. The max/min range values are: median +/- N * sigmaMad"
67 ", where N is this config value.",
68 default=3.0,
69 )
70 colorMapType = ChoiceField[str](
71 doc="Type of color map to use for the color bar. Options: sequential, divergent, userDefined.",
72 allowed={cmType: cmType for cmType in ("sequential", "divergent")},
73 default="divergent",
74 )
75 colorMap = ListField[str](
76 doc="List of hexidecimal colors for a user-defined color map.",
77 optional=True,
78 )
79 showOutliers = Field[bool](
80 doc="Show the outliers on the plot. "
81 "Outliers are values whose absolute value is > colorBarRange * sigmaMAD.",
82 default=True,
83 )
84 showNaNs = Field[bool](doc="Show the NaNs on the plot.", default=True)
85 labelTracts = Field[bool](doc="Label the tracts.", default=False)
87 def getInputSchema(self, **kwargs) -> KeyedDataSchema:
88 base = []
89 base.append(("z", Vector))
90 base.append(("tract", Vector))
91 return base
93 def __call__(self, data: KeyedData, **kwargs) -> Mapping[str, Figure] | Figure:
94 self._validateInput(data, **kwargs)
95 return self.makePlot(data, **kwargs)
97 def _validateInput(self, data: KeyedData, **kwargs) -> None:
98 """NOTE currently can only check that something is not a Scalar, not
99 check that the data is consistent with Vector
100 """
101 needed = self.getInputSchema(**kwargs)
102 if remainder := {key.format(**kwargs) for key, _ in needed} - {
103 key.format(**kwargs) for key in data.keys()
104 }:
105 raise ValueError(f"Task needs keys {remainder} but they were not found in input")
106 for name, typ in needed:
107 isScalar = issubclass((colType := type(data[name.format(**kwargs)])), Scalar)
108 if isScalar and typ != Scalar:
109 raise ValueError(f"Data keyed by {name} has type {colType} but action requires type {typ}")
111 def _getAxesLimits(self, xs: list, ys: list) -> tuple(list, list):
112 """Get the x and y axes limits in degrees.
114 Parameters
115 ----------
116 xs : `list`
117 X coordinates for the tracts to plot.
118 ys : `list`
119 Y coordinates for the tracts to plot.
121 Returns
122 -------
123 xlim : `list`
124 Minimun and maximum x axis values.
125 ylim : `list`
126 Minimun and maximum y axis values.
127 """
129 # Add some blank space on the edges of the plot.
130 xlim = [np.nanmin(xs) - 5, np.nanmax(xs) + 5]
131 ylim = [np.nanmin(ys) - 5, np.nanmax(ys) + 5]
133 # Limit to only show real RA/Dec values.
134 if xlim[0] < 0.0:
135 xlim[0] = 0.0
136 if xlim[1] > 360.0:
137 xlim[1] = 360.0
138 if ylim[0] < -90.0:
139 ylim[0] = -90.0
140 if ylim[1] > 90.0:
141 ylim[1] = 90.0
143 return (xlim, ylim)
145 def _getMaxOutlierVals(self, multiplier: float, tracts: list, values: list, outlierInds: list) -> str:
146 """Get the 5 largest outlier values in a string.
148 Parameters
149 ----------
150 multiplier : `float`
151 Select values whose absolute value is > multiplier * sigmaMAD.
152 tracts : `list`
153 All the tracts.
154 values : `list`
155 All the metric values.
156 outlierInds : `list`
157 Indicies of outlier values.
159 Returns
160 -------
161 text : `str`
162 A string containing the 5 tracts with the largest outlier values.
163 """
164 text = f"Tracts with |value| > {multiplier}" + r"$\sigma_{MAD}$" + ": "
165 if len(outlierInds) > 0:
166 outlierValues = np.array(values)[outlierInds]
167 outlierTracts = np.array(tracts)[outlierInds]
168 # Sort values in descending (-) absolute value order discounting
169 # NaNs.
170 maxInds = np.argsort(-np.abs(outlierValues))
171 # Show up to five values on the plot.
172 for ind in maxInds[:5]:
173 val = outlierValues[ind]
174 tract = outlierTracts[ind]
175 text += f"{tract}, {val:.3}; "
176 # Remove the final trailing comma and whitespace.
177 text = text[:-2]
178 else:
179 text += "None"
181 return text
183 def makePlot(
184 self,
185 data: KeyedData,
186 plotInfo: Optional[Mapping[str, str]] = None,
187 **kwargs,
188 ) -> Figure:
189 """Make a WholeSkyPlot of the given data.
191 Parameters
192 ----------
193 data : `KeyedData`
194 The catalog to plot the points from.
195 plotInfo : `dict`
196 A dictionary of information about the data being plotted with keys:
198 ``"run"``
199 The output run for the plots (`str`).
200 ``"skymap"``
201 The type of skymap used for the data (`str`).
202 ``"filter"``
203 The filter used for this data (`str`).
204 ``"tract"``
205 The tract that the data comes from (`str`).
207 Returns
208 -------
209 `pipeBase.Struct` containing:
210 skyPlot : `matplotlib.figure.Figure`
211 The resulting figure.
214 Examples
215 --------
216 An example of the plot produced from this code is here:
218 .. image:: /_static/analysis_tools/wholeSkyPlotExample.png
220 For a detailed example of how to make a plot from the command line
221 please see the
222 :ref:`getting started guide<analysis-tools-getting-started>`.
223 """
224 skymap = kwargs["skymap"]
225 if plotInfo is None:
226 plotInfo = {}
228 # Prevent Bands in the plot info showing a list of bands.
229 # If bands is a list, it implies that parameterizedBand=False,
230 # and that the metric is not band-specific.
231 if "bands" in plotInfo:
232 if isinstance(plotInfo["bands"], list):
233 plotInfo["bands"] = "N/A"
235 colorMap = self.colorMap
236 match self.colorMapType:
237 case "sequential":
238 if colorMap is None:
239 colorMap = mkColormap(["#F5F5F5", "#5AB4AC", "#284D48"])
240 norm = None
241 case "divergent":
242 if colorMap is None:
243 colorMap = mkColormap(["#9A6E3A", "#C6A267", "#A9A9A9", "#4F938B", "#2C665A"])
244 norm = CenteredNorm()
246 # Create patches using the corners of each tract.
247 patches = []
248 colBarVals = []
249 tracts = []
250 ras = []
251 decs = []
252 mid_ras = []
253 mid_decs = []
254 for i, tract in enumerate(data["tract"]):
255 corners = getTractCorners(skymap, tract)
256 patches.append(Polygon(corners, closed=True))
257 colBarVals.append(data["z"][i])
258 tracts.append(tract)
259 ras.append(corners[0][0])
260 decs.append(corners[0][1])
261 mid_ras.append((corners[0][0] + corners[1][0]) / 2)
262 mid_decs.append((corners[0][1] + corners[2][1]) / 2)
264 # Setup figure.
265 fig, ax = plt.subplots(1, 1, figsize=self.figureSize, dpi=self.dpi)
266 if self.autoAxesLimits:
267 xlim, ylim = self._getAxesLimits(ras, decs)
268 else:
269 xlim, ylim = self.xLimits, self.yLimits
270 ax.set_xlim(xlim)
271 ax.set_ylim(ylim)
272 ax.set_xlabel(self.xAxisLabel)
273 ax.set_ylabel(self.yAxisLabel)
274 ax.invert_xaxis()
276 # Add colored patches showing tract metric values.
277 patchCollection = PatchCollection(patches, cmap=colorMap, norm=norm)
278 ax.add_collection(patchCollection)
280 # Define color bar range.
281 med = np.nanmedian(colBarVals)
282 sigmaMad = nanSigmaMad(colBarVals)
283 if self.colorBarMin is not None:
284 vmin = np.float64(self.colorBarMin)
285 else:
286 vmin = med - self.colorBarRange * sigmaMad
287 if self.colorBarMax is not None:
288 vmax = np.float64(self.colorBarMax)
289 else:
290 vmax = med + self.colorBarRange * sigmaMad
292 # Note tracts with metrics outside (vmin, vmax) as outliers.
293 outlierInds = np.where((colBarVals < vmin) | (colBarVals > vmax))[0]
295 # Initialize legend handles.
296 handles = []
298 if self.showOutliers:
299 # Plot the outlier patches.
300 outlierPatches = []
301 if len(outlierInds) > 0:
302 for ind in outlierInds:
303 outlierPatches.append(patches[ind])
304 outlierPatchCollection = PatchCollection(
305 outlierPatches,
306 cmap=colorMap,
307 norm=norm,
308 facecolors="none",
309 edgecolors="k",
310 linewidths=0.5,
311 zorder=100,
312 )
313 ax.add_collection(outlierPatchCollection)
314 # Add legend information.
315 outlierPatch = Patch(
316 facecolor="none",
317 edgecolor="k",
318 linewidth=0.5,
319 label="Outlier",
320 )
321 handles.append(outlierPatch)
323 if self.showNaNs:
324 # Plot tracts with NaN metric values.
325 nanInds = np.where(~np.isfinite(colBarVals))[0]
326 nanPatches = []
327 if len(nanInds) > 0:
328 for ind in nanInds:
329 nanPatches.append(patches[ind])
330 nanPatchCollection = PatchCollection(
331 nanPatches,
332 cmap=None,
333 norm=norm,
334 facecolors="white",
335 edgecolors="grey",
336 linestyles="dotted",
337 linewidths=0.5,
338 zorder=100,
339 )
340 ax.add_collection(nanPatchCollection)
341 # Add legend information.
342 nanPatch = Patch(
343 facecolor="white",
344 edgecolor="grey",
345 linestyle="dotted",
346 linewidth=0.5,
347 label="NaN",
348 )
349 handles.append(nanPatch)
351 if len(handles) > 0:
352 fig.legend(handles=handles)
354 if self.labelTracts:
355 # Label the tracts
356 for i, tract in enumerate(tracts):
357 ax.text(
358 mid_ras[i],
359 mid_decs[i],
360 f"{tract}",
361 ha="center",
362 va="center",
363 fontsize=2,
364 alpha=0.7,
365 zorder=100,
366 )
368 if self.showOutliers:
369 # Add text boxes to show the number of tracts, number of NaNs,
370 # median, sigma MAD, and the five largest outlier values.
371 outlierText = self._getMaxOutlierVals(self.colorBarRange, tracts, colBarVals, outlierInds)
373 # Make vertical text spacing readable for different figure sizes.
374 multiplier = 3.5 / self.figureSize[1]
375 verticalSpacing = 0.028 * multiplier
376 fig.text(
377 0.01,
378 0.01 + 3 * verticalSpacing,
379 f"Num tracts: {len(tracts)}",
380 transform=fig.transFigure,
381 fontsize=8,
382 alpha=0.7,
383 )
384 if self.showNaNs:
385 fig.text(
386 0.01,
387 0.01 + 2 * verticalSpacing,
388 f"Num nans: {len(nanInds)}",
389 transform=fig.transFigure,
390 fontsize=8,
391 alpha=0.7,
392 )
393 fig.text(
394 0.01,
395 0.01 + verticalSpacing,
396 f"Median: {med:.3f}; " + r"$\sigma_{MAD}$" + f": {sigmaMad:.3f}",
397 transform=fig.transFigure,
398 fontsize=8,
399 alpha=0.7,
400 )
401 if self.showOutliers:
402 fig.text(0.01, 0.01, outlierText, transform=fig.transFigure, fontsize=8, alpha=0.7)
404 # Truncate the color range to (vmin, vmax).
405 colorBarVals = np.clip(np.array(colBarVals), vmin, vmax)
406 patchCollection.set_array(colorBarVals)
407 # Make the color bar with a metric label.
408 cbar = plt.colorbar(
409 patchCollection,
410 ax=ax,
411 shrink=0.7,
412 extend="both",
413 location="top",
414 orientation="horizontal",
415 )
416 cbarText = self.zAxisLabel.format_map(kwargs)
417 if "zUnit" in data and data["zUnit"] != "":
418 cbarText += f" ({data['zUnit']})"
419 text = cbar.ax.text(
420 0.5,
421 0.5,
422 cbarText,
423 transform=cbar.ax.transAxes,
424 ha="center",
425 va="center",
426 fontsize=10,
427 zorder=100,
428 )
429 text.set_path_effects([pathEffects.Stroke(linewidth=3, foreground="w"), pathEffects.Normal()])
431 # Finalize plot appearance.
432 ax.grid()
433 ax.set_axisbelow(True)
434 ax.set_aspect("equal")
435 fig = plt.gcf()
436 fig = addPlotInfo(fig, plotInfo)
437 plt.subplots_adjust(left=0.08, right=0.97, top=0.8, bottom=0.17, wspace=0.35)
439 return fig