Coverage for python / lsst / analysis / tools / actions / plot / wholeSkyPlot.py: 14%
241 statements
« prev ^ index » next coverage.py v7.13.2, created at 2026-02-14 02:49 -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",)
26import importlib.resources as importResources
27import json
28from typing import Mapping, Optional
30import lsst.analysis.tools
31import matplotlib.patheffects as pathEffects
32import numpy as np
33import yaml
34from lsst.pex.config import ChoiceField, Field, ListField
35from lsst.utils.plotting import (
36 accent_color,
37 divergent_cmap,
38 make_figure,
39 set_rubin_plotstyle,
40 stars_cmap,
41 stars_color,
42)
43from matplotlib import gridspec
44from matplotlib.collections import PatchCollection
45from matplotlib.colors import CenteredNorm
46from matplotlib.figure import Figure
47from matplotlib.patches import Patch, Polygon
49from ...interfaces import KeyedData, KeyedDataSchema, PlotAction, Scalar, Vector
50from ...math import nanSigmaMad
51from ...utils import getTractCorners
52from .plotUtils import addPlotInfo
55class AnnotatedFigure(Figure):
56 metadata: dict
59class WholeSkyPlot(PlotAction):
60 """Plots the on sky distribution of a parameter.
62 Plots the values of the parameter given for the z axis
63 according to the positions given for x and y. Optimised
64 for use with RA and Dec. Also calculates some basic
65 statistics and includes those on the plot.
67 The default axes limits and figure size were chosen to plot HSC PDR2.
68 """
70 xAxisLabel = Field[str](doc="Label to use for the x axis.", default="RA (degrees)")
71 yAxisLabel = Field[str](doc="Label to use for the y axis.", default="Dec (degrees)")
72 zAxisLabel = Field[str](doc="Label to use for the z axis.", default="")
73 autoAxesLimits = Field[bool](doc="Find axes limits automatically.", default=True)
74 xLimits = ListField[float](doc="Plotting limits for the x axis.", default=[-5.0, 365.0])
75 yLimits = ListField[float](doc="Plotting limits for the y axis.", default=[-10.0, 60.0])
76 autoAxesLimits = Field[bool](doc="Find axes limits automatically.", default=True)
77 colorBarMin = Field[float](doc="The minimum value of the color bar.", optional=True)
78 colorBarMax = Field[float](doc="The minimum value of the color bar.", optional=True)
79 colorBarRange = Field[float](
80 doc="The multiplier for the color bar range. The max/min range values are: median +/- N * sigmaMad"
81 ", where N is this config value.",
82 default=3.0,
83 )
84 colorMapType = ChoiceField[str](
85 doc="Type of color map to use for the color bar. Options: sequential, divergent, userDefined.",
86 allowed={cmType: cmType for cmType in ("sequential", "divergent")},
87 default="divergent",
88 )
89 colorMap = ListField[str](
90 doc="List of hexidecimal colors for a user-defined color map.",
91 optional=True,
92 )
93 showOutliers = Field[bool](
94 doc="Show the outliers on the plot. "
95 "Outliers are values whose absolute value is > colorBarRange * sigmaMAD.",
96 default=True,
97 )
98 showNaNs = Field[bool](doc="Show the NaNs on the plot.", default=True)
99 labelTracts = Field[bool](doc="Label the tracts.", default=False)
101 addThresholds = Field[bool](
102 doc="Read in the predefined thresholds and indicate them on the histogram.",
103 default=True,
104 )
106 def getInputSchema(self, **kwargs) -> KeyedDataSchema:
107 base = []
108 base.append(("z", Vector))
109 base.append(("tract", Vector))
110 return base
112 def __call__(self, data: KeyedData, **kwargs) -> Mapping[str, Figure] | Figure:
113 self._validateInput(data, **kwargs)
114 return self.makePlot(data, **kwargs)
116 def _validateInput(self, data: KeyedData, **kwargs) -> None:
117 """NOTE currently can only check that something is not a Scalar, not
118 check that the data is consistent with Vector
119 """
120 needed = self.getInputSchema(**kwargs)
121 if remainder := {key.format(**kwargs) for key, _ in needed} - {
122 key.format(**kwargs) for key in data.keys()
123 }:
124 raise ValueError(f"Task needs keys {remainder} but they were not found in input")
125 for name, typ in needed:
126 isScalar = issubclass((colType := type(data[name.format(**kwargs)])), Scalar)
127 if isScalar and typ != Scalar:
128 raise ValueError(f"Data keyed by {name} has type {colType} but action requires type {typ}")
130 def _getAxesLimits(self, xs: list, ys: list) -> tuple(list, list):
131 """Get the x and y axes limits in degrees.
133 Parameters
134 ----------
135 xs : `list`
136 X coordinates for the tracts to plot.
137 ys : `list`
138 Y coordinates for the tracts to plot.
140 Returns
141 -------
142 xlim : `list`
143 Minimun and maximum x axis values.
144 ylim : `list`
145 Minimun and maximum y axis values.
146 """
148 # Add some blank space on the edges of the plot.
149 xlim = [np.nanmin(xs) - 5, np.nanmax(xs) + 5]
150 ylim = [np.nanmin(ys) - 5, np.nanmax(ys) + 5]
152 # Limit to only show real RA/Dec values.
153 if xlim[0] < 0.0:
154 xlim[0] = 0.0
155 if xlim[1] > 360.0:
156 xlim[1] = 360.0
157 if ylim[0] < -90.0:
158 ylim[0] = -90.0
159 if ylim[1] > 90.0:
160 ylim[1] = 90.0
162 return (xlim, ylim)
164 def _getMaxOutlierVals(self, multiplier: float, tracts: list, values: list, outlierInds: list) -> str:
165 """Get the 5 largest outlier values in a string.
167 Parameters
168 ----------
169 multiplier : `float`
170 Select values whose absolute value is > multiplier * sigmaMAD.
171 tracts : `list`
172 All the tracts.
173 values : `list`
174 All the metric values.
175 outlierInds : `list`
176 Indicies of outlier values.
178 Returns
179 -------
180 text : `str`
181 A string containing the 10 tracts with the largest outlier values.
182 """
183 if self.addThresholds:
184 text = "Tracts with value outside thresholds: "
185 else:
186 text = f"Tracts with |value| > {multiplier}" + r"$\sigma_{MAD}$" + ": "
187 if len(outlierInds) > 0:
188 outlierValues = np.array(values)[outlierInds]
189 outlierTracts = np.array(tracts)[outlierInds]
190 # Sort values in descending (-) absolute value order discounting
191 # NaNs.
192 maxInds = np.argsort(-np.abs(outlierValues))
193 # Show up to ten values on the plot.
194 for ind in maxInds[:10]:
195 val = outlierValues[ind]
196 tract = outlierTracts[ind]
197 text += f"{tract}, {val:.3}; "
198 # Remove the final trailing comma and whitespace.
199 text = text[:-2]
200 else:
201 text += "None"
203 return text
205 def makePlot(
206 self,
207 data: KeyedData,
208 plotInfo: Optional[Mapping[str, str]] = None,
209 **kwargs,
210 ) -> AnnotatedFigure:
211 """Make a WholeSkyPlot of the given data.
213 Parameters
214 ----------
215 data : `KeyedData`
216 The catalog to plot the points from.
217 plotInfo : `dict`
218 A dictionary of information about the data being plotted with keys:
220 ``"run"``
221 The output run for the plots (`str`).
222 ``"skymap"``
223 The type of skymap used for the data (`str`).
224 ``"filter"``
225 The filter used for this data (`str`).
226 ``"tract"``
227 The tract that the data comes from (`str`).
229 Returns
230 -------
231 `pipeBase.Struct` containing:
232 skyPlot : `matplotlib.figure.Figure`
233 The resulting figure.
236 Examples
237 --------
238 An example of the plot produced from this code is here:
240 .. image:: /_static/analysis_tools/wholeSkyPlotExample.png
242 For a detailed example of how to make a plot from the command line
243 please see the
244 :ref:`getting started guide<analysis-tools-getting-started>`.
245 """
246 skymap = kwargs["skymap"]
247 if plotInfo is None:
248 plotInfo = {}
250 if self.addThresholds:
251 metricThresholdFile = importResources.read_text(lsst.analysis.tools, "metricInformation.yaml")
252 metricDefs = yaml.safe_load(metricThresholdFile)
254 # Prevent Bands in the plot info showing a list of bands.
255 # If bands is a list, it implies that parameterizedBand=False,
256 # and that the metric is not band-specific.
257 if "bands" in plotInfo:
258 if isinstance(plotInfo["bands"], list):
259 plotInfo["bands"] = "N/A"
261 colorMap = self.colorMap
262 match self.colorMapType:
263 case "sequential":
264 if colorMap is None:
265 colorMap = stars_cmap()
266 outlierColor = "red"
267 norm = None
268 case "divergent":
269 if colorMap is None:
270 colorMap = divergent_cmap()
271 outlierColor = "fuchsia"
272 norm = CenteredNorm()
274 # Create patches using the corners of each tract.
275 patches = []
276 colBarVals = []
277 tracts = []
278 ras = []
279 decs = []
280 mid_ras = []
281 mid_decs = []
282 for i, tract in enumerate(data["tract"]):
283 corners = getTractCorners(skymap, tract)
284 patches.append(Polygon(corners, closed=True))
285 colBarVals.append(data["z"][i])
286 tracts.append(tract)
287 ras.append(corners[0][0])
288 decs.append(corners[0][1])
289 mid_ras.append((corners[0][0] + corners[1][0]) / 2)
290 mid_decs.append((corners[0][1] + corners[2][1]) / 2)
292 # Setup figure.
293 fig: AnnotatedFigure = make_figure(dpi=300, figsize=(12, 3.5))
294 set_rubin_plotstyle()
295 gs = gridspec.GridSpec(1, 4)
296 ax = fig.add_subplot(gs[:3])
297 # Add colored patches showing tract metric values.
298 patchCollection = PatchCollection(patches, cmap=colorMap, norm=norm)
299 ax.add_collection(patchCollection)
301 # Define color bar range.
302 if np.sum(np.isfinite(colBarVals)) > 0:
303 med = np.nanmedian(colBarVals)
304 else:
305 med = np.nan
306 sigmaMad = nanSigmaMad(colBarVals)
307 if self.colorBarMin is not None:
308 vmin = np.float64(self.colorBarMin)
309 else:
310 vmin = med - self.colorBarRange * sigmaMad
311 if self.colorBarMax is not None:
312 vmax = np.float64(self.colorBarMax)
313 else:
314 vmax = med + self.colorBarRange * sigmaMad
316 dataName = self.zAxisLabel.format_map(kwargs)
317 colBarVals = np.array(colBarVals)
318 if self.addThresholds and dataName in metricDefs:
319 if "lowThreshold" in metricDefs[dataName].keys():
320 lowThreshold = metricDefs[dataName]["lowThreshold"]
321 else:
322 lowThreshold = np.nan
323 if "highThreshold" in metricDefs[dataName].keys():
324 highThreshold = metricDefs[dataName]["highThreshold"]
325 else:
326 highThreshold = np.nan
327 outlierInds = np.where((colBarVals < lowThreshold) | (colBarVals > highThreshold))[0]
328 else:
329 # Note tracts with metrics outside (vmin, vmax) as outliers.
330 outlierInds = np.where((colBarVals < vmin) | (colBarVals > vmax))[0]
332 # Initialize legend handles.
333 handles = []
335 if self.showOutliers:
336 # Plot the outlier patches.
337 outlierPatches = []
338 if len(outlierInds) > 0:
339 for ind in outlierInds:
340 outlierPatches.append(patches[ind])
341 outlierPatchCollection = PatchCollection(
342 outlierPatches,
343 cmap=colorMap,
344 norm=norm,
345 facecolors="none",
346 edgecolors=outlierColor,
347 linewidths=0.5,
348 zorder=100,
349 )
350 ax.add_collection(outlierPatchCollection)
351 # Add legend information.
352 outlierPatch = Patch(
353 facecolor="none",
354 edgecolor=outlierColor,
355 linewidth=0.5,
356 label="Outlier",
357 )
358 handles.append(outlierPatch)
360 if self.showNaNs:
361 # Plot tracts with NaN metric values.
362 nanInds = np.where(~np.isfinite(colBarVals))[0]
363 nanPatches = []
364 if len(nanInds) > 0:
365 for ind in nanInds:
366 nanPatches.append(patches[ind])
367 nanPatchCollection = PatchCollection(
368 nanPatches,
369 cmap=None,
370 norm=norm,
371 facecolors="white",
372 edgecolors="grey",
373 linestyles="dotted",
374 linewidths=0.5,
375 zorder=100,
376 )
377 ax.add_collection(nanPatchCollection)
378 # Add legend information.
379 nanPatch = Patch(
380 facecolor="white",
381 edgecolor="grey",
382 linestyle="dotted",
383 linewidth=0.5,
384 label="NaN",
385 )
386 handles.append(nanPatch)
388 if len(handles) > 0:
389 fig.legend(handles=handles)
391 if self.labelTracts:
392 # Label the tracts
393 for i, tract in enumerate(tracts):
394 ax.text(
395 mid_ras[i],
396 mid_decs[i],
397 f"{tract}",
398 ha="center",
399 va="center",
400 fontsize=2,
401 alpha=0.7,
402 zorder=100,
403 )
405 ax.set_aspect("equal")
406 axPos = ax.get_position()
407 ax1 = fig.add_axes([0.73, 0.25, 0.20, 0.47])
409 if np.sum(np.isfinite(data["z"])) > 0:
410 ax1.hist(data["z"], bins=len(data["z"] / 10), color=stars_color(), histtype="step")
411 else:
412 ax1.text(0.5, 0.5, "Data all NaN/Inf")
413 ax1.set_xlabel("Metric Values")
414 ax1.set_ylabel("Number")
415 ax1.yaxis.set_label_position("right")
416 ax1.yaxis.tick_right()
418 if self.addThresholds and dataName in metricDefs:
419 # Check the thresholds are finite and set them to
420 # the min/max of the data if they aren't to calculate
421 # the x range of the plot
422 if np.isfinite(lowThreshold):
423 ax1.axvline(lowThreshold, color=accent_color())
424 else:
425 lowThreshold = np.nanmin(colBarVals)
426 if np.isfinite(highThreshold):
427 ax1.axvline(highThreshold, color=accent_color())
428 else:
429 highThreshold = np.nanmax(colBarVals)
431 widthThreshold = highThreshold - lowThreshold
432 upperLim = highThreshold + 0.5 * widthThreshold
433 lowerLim = lowThreshold - 0.5 * widthThreshold
434 ax1.set_xlim(lowerLim, upperLim)
435 numOutside = np.sum(((data["z"] > upperLim) | (data["z"] < lowerLim)))
436 ax1.set_title("Outside plot limits: " + str(numOutside))
438 else:
439 if vmin != vmax and np.isfinite(vmin) and np.isfinite(vmax):
440 ax1.set_xlim(vmin, vmax)
442 if self.autoAxesLimits:
443 xlim, ylim = self._getAxesLimits(ras, decs)
444 else:
445 xlim, ylim = self.xLimits, self.yLimits
446 ax.set_xlim(xlim)
447 ax.set_ylim(ylim)
448 ax.set_xlabel(self.xAxisLabel)
449 ax.set_ylabel(self.yAxisLabel)
450 ax.invert_xaxis()
452 if self.showOutliers:
453 # Add text boxes to show the number of tracts, number of NaNs,
454 # median, sigma MAD, and the five largest outlier values.
455 outlierText = self._getMaxOutlierVals(self.colorBarRange, tracts, colBarVals, outlierInds)
456 # Make vertical text spacing readable for different figure sizes.
457 multiplier = 3.5 / fig.get_size_inches()[1]
458 verticalSpacing = 0.028 * multiplier
459 fig.text(
460 0.01,
461 0.01 + 3 * verticalSpacing,
462 f"Num tracts: {len(tracts)}",
463 transform=fig.transFigure,
464 fontsize=8,
465 alpha=0.7,
466 )
467 if self.showNaNs:
468 fig.text(
469 0.01,
470 0.01 + 2 * verticalSpacing,
471 f"Num nans: {len(nanInds)}",
472 transform=fig.transFigure,
473 fontsize=8,
474 alpha=0.7,
475 )
476 fig.text(
477 0.01,
478 0.01 + verticalSpacing,
479 f"Median: {med:.3f}; " + r"$\sigma_{MAD}$" + f": {sigmaMad:.3f}",
480 transform=fig.transFigure,
481 fontsize=8,
482 alpha=0.7,
483 )
484 if self.showOutliers:
485 fig.text(0.01, 0.01, outlierText, transform=fig.transFigure, fontsize=8, alpha=0.7)
487 # Truncate the color range to (vmin, vmax).
488 if vmin != vmax and np.isfinite(vmin) and np.isfinite(vmax):
489 colBarVals = np.clip(np.array(colBarVals), vmin, vmax)
490 patchCollection.set_array(colBarVals)
491 # Make the color bar with a metric label.
492 axPos = ax.get_position()
493 cax = fig.add_axes([0.084, axPos.y1 + 0.02, 0.62, 0.07])
494 fig.colorbar(
495 patchCollection,
496 cax=cax,
497 shrink=0.7,
498 extend="both",
499 location="top",
500 orientation="horizontal",
501 )
502 cbarText = "Metric Values"
504 text = cax.text(
505 0.5,
506 0.5,
507 cbarText,
508 transform=cax.transAxes,
509 ha="center",
510 va="center",
511 fontsize=10,
512 zorder=100,
513 )
514 text.set_path_effects([pathEffects.Stroke(linewidth=3, foreground="w"), pathEffects.Normal()])
516 # Finalize plot appearance.
517 ax.grid()
518 ax.set_axisbelow(True)
519 addPlotInfo(fig, plotInfo)
520 fig.subplots_adjust(left=0.08, right=0.92, top=0.8, bottom=0.17, wspace=0.05)
521 titleText = self.zAxisLabel.format_map(kwargs)
522 if "zUnit" in data and data["zUnit"] != "":
523 titleText += f" ({data['zUnit']})"
524 fig.suptitle("Metric: " + titleText, fontsize=20)
526 # This saves metadata in the PNG that allows the plot-navigator
527 # to provide tract numbers and metric values on mouseover.
528 #
529 # PNG metadata is a set of string keys and string values.
530 # The WholeSkyPlot stores two keys:
531 # - label: the string describing the regions ('tract')
532 # - boxes, JSON string of a list of per-region dictionaries,
533 # where each dictionary has fields:
534 # - min_x, max_x, min_y, max_y for the pixel coordinates of
535 # the four corners of the region
536 # - id: the identifier of the region (e.g. tract number)
537 # - value: the region's metric, as a string.
538 #
539 def make_patch_md(patch, id_field, value, ax):
540 path = ax.transData.transform_path(patch.get_path())
541 x_path = [int(x) for x in path.vertices[:, 0].tolist()]
542 y_path = [int(y) for y in path.vertices[:, 1].tolist()]
543 return {
544 "min_x": min(x_path),
545 "max_x": max(x_path),
546 "min_y": min(y_path),
547 "max_y": max(y_path),
548 "id": f"{id_field}",
549 "value": f"{value:.3}",
550 }
552 # After ax.set_aspect(), the figure needs to be drawn for the axes
553 # transformations to be updated to the right values.
554 fig.canvas.draw_idle()
556 patch_coordinate_entries = [
557 make_patch_md(patch, tract, value, ax)
558 for (patch, tract, value) in zip(patches, tracts, colBarVals)
559 ]
561 fig.metadata = {"label": "Tract", "boxes": json.dumps(patch_coordinate_entries)}
563 return fig