Coverage for python/lsst/analysis/tools/actions/plot/wholeSkyPlot.py: 13%
265 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-08 00:04 +0000
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 copy
27import importlib.resources as importResources
28import json
29from collections.abc import Mapping
31import matplotlib.patheffects as pathEffects
32import numpy as np
33import skyproj
34import yaml
35from matplotlib import colorizer, gridspec
36from matplotlib.colors import Normalize
37from matplotlib.figure import Figure
38from matplotlib.patches import Rectangle
40import lsst.analysis.tools
41from lsst.pex.config import ChoiceField, Field, ListField
42from lsst.utils.plotting import (
43 accent_color,
44 divergent_cmap,
45 make_figure,
46 set_rubin_plotstyle,
47 stars_cmap,
48 stars_color,
49)
51from ...interfaces import KeyedData, KeyedDataSchema, PlotAction, Scalar, Vector
52from ...math import nanSigmaMad
53from ...utils import getTractCorners
54from .plotUtils import addPlotInfo
57class AnnotatedFigure(Figure):
58 metadata: dict
61class WholeSkyPlot(PlotAction):
62 """Plots the on sky distribution of a parameter.
64 Plots the values of the parameter given for the z axis
65 according to the positions given for x and y. Optimised
66 for use with RA and Dec. Also calculates some basic
67 statistics and includes those on the plot.
69 The default axes limits and figure size were chosen to plot HSC PDR2.
70 """
72 xAxisLabel = Field[str](doc="Label to use for the x axis.", default="R.A.")
73 yAxisLabel = Field[str](doc="Label to use for the y axis.", default="Dec.")
74 zAxisLabel = Field[str](doc="Label to use for the z axis.", default="")
75 autoAxesLimits = Field[bool](doc="Find axes limits automatically.", default=True)
76 xLimits = ListField[float](doc="Plotting limits for the x axis.", default=[-5.0, 365.0])
77 yLimits = ListField[float](doc="Plotting limits for the y axis.", default=[-10.0, 60.0])
78 colorBarMin = Field[float](doc="The minimum value of the color bar.", optional=True)
79 colorBarMax = Field[float](doc="The minimum value of the color bar.", optional=True)
80 colorBarRange = Field[float](
81 doc="The multiplier for the color bar range. The max/min range values are: median +/- N * sigmaMad"
82 ", where N is this config value.",
83 default=3.0,
84 )
85 colorMapType = ChoiceField[str](
86 doc="Type of color map to use for the color bar. Options: sequential, divergent, userDefined.",
87 allowed={cmType: cmType for cmType in ("sequential", "divergent")},
88 default="divergent",
89 )
90 colorMap = ListField[str](
91 doc="List of hexidecimal colors for a user-defined color map.",
92 optional=True,
93 )
94 showOutliers = Field[bool](
95 doc="Show the outliers on the plot. "
96 "Outliers are values whose absolute value is > colorBarRange * sigmaMAD.",
97 default=True,
98 )
99 showNaNs = Field[bool](doc="Show the NaNs on the plot.", default=True)
100 labelTracts = Field[bool](doc="Label the tracts.", default=False)
102 addThresholds = Field[bool](
103 doc="Read in the predefined thresholds and indicate them on the histogram.",
104 default=True,
105 )
107 publicationStyle = Field[bool](
108 doc="Make a simplified plot for publication use.",
109 default=False,
110 )
112 fixAroundZero = Field[bool](
113 doc="Fix the colorbar to be symmetric around zero.",
114 default=False,
115 )
117 def getInputSchema(self, **kwargs) -> KeyedDataSchema:
118 base = []
119 base.append(("z", Vector))
120 base.append(("tract", Vector))
121 return base
123 def __call__(self, data: KeyedData, **kwargs) -> Mapping[str, Figure] | Figure:
124 self._validateInput(data, **kwargs)
125 return self.makePlot(data, **kwargs)
127 def _validateInput(self, data: KeyedData, **kwargs) -> None:
128 """NOTE currently can only check that something is not a Scalar, not
129 check that the data is consistent with Vector
130 """
131 needed = self.getInputSchema(**kwargs)
132 if remainder := {key.format(**kwargs) for key, _ in needed} - {
133 key.format(**kwargs) for key in data.keys()
134 }:
135 raise ValueError(f"Task needs keys {remainder} but they were not found in input")
136 for name, typ in needed:
137 isScalar = issubclass((colType := type(data[name.format(**kwargs)])), Scalar)
138 if isScalar and typ != Scalar:
139 raise ValueError(f"Data keyed by {name} has type {colType} but action requires type {typ}")
141 def _getAxesLimits(self, xs: list, ys: list) -> tuple(list, list):
142 """Get the x and y axes limits in degrees.
144 Parameters
145 ----------
146 xs : `list`
147 X coordinates for the tracts to plot.
148 ys : `list`
149 Y coordinates for the tracts to plot.
151 Returns
152 -------
153 xlim : `list`
154 Minimun and maximum x axis values.
155 ylim : `list`
156 Minimun and maximum y axis values.
157 """
159 # Add some blank space on the edges of the plot.
160 xlim = [np.nanmin(xs) - 5, np.nanmax(xs) + 5]
161 ylim = [np.nanmin(ys) - 5, np.nanmax(ys) + 5]
163 # Limit to only show real RA/Dec values.
164 if xlim[0] < 0.0:
165 xlim[0] = 0.0
166 if xlim[1] > 360.0:
167 xlim[1] = 360.0
168 if ylim[0] < -90.0:
169 ylim[0] = -90.0
170 if ylim[1] > 90.0:
171 ylim[1] = 90.0
173 return (xlim, ylim)
175 def _getMaxOutlierVals(self, multiplier: float, tracts: list, values: list, outlierInds: list) -> str:
176 """Get the 5 largest outlier values in a string.
178 Parameters
179 ----------
180 multiplier : `float`
181 Select values whose absolute value is > multiplier * sigmaMAD.
182 tracts : `list`
183 All the tracts.
184 values : `list`
185 All the metric values.
186 outlierInds : `list`
187 Indicies of outlier values.
189 Returns
190 -------
191 text : `str`
192 A string containing the 10 tracts with the largest outlier values.
193 """
194 if self.addThresholds:
195 text = "Tracts with value outside thresholds: "
196 else:
197 text = f"Tracts with |value| > {multiplier}" + r"$\sigma_{MAD}$" + ": "
198 if len(outlierInds) > 0:
199 outlierValues = np.array(values)[outlierInds]
200 outlierTracts = np.array(tracts)[outlierInds]
201 # Sort values in descending (-) absolute value order discounting
202 # NaNs.
203 maxInds = np.argsort(-np.abs(outlierValues))
204 # Show up to ten values on the plot.
205 for ind in maxInds[:10]:
206 val = outlierValues[ind]
207 tract = outlierTracts[ind]
208 text += f"{tract}, {val:.3}; "
209 # Remove the final trailing comma and whitespace.
210 text = text[:-2]
211 else:
212 text += "None"
214 return text
216 def makePlot(
217 self,
218 data: KeyedData,
219 plotInfo: Mapping[str, str] | None = None,
220 **kwargs,
221 ) -> AnnotatedFigure:
222 """Make a WholeSkyPlot of the given data.
224 Parameters
225 ----------
226 data : `KeyedData`
227 The catalog to plot the points from.
228 plotInfo : `dict`
229 A dictionary of information about the data being plotted with keys:
231 ``"run"``
232 The output run for the plots (`str`).
233 ``"skymap"``
234 The type of skymap used for the data (`str`).
235 ``"filter"``
236 The filter used for this data (`str`).
237 ``"tract"``
238 The tract that the data comes from (`str`).
240 Returns
241 -------
242 `pipeBase.Struct` containing:
243 skyPlot : `matplotlib.figure.Figure`
244 The resulting figure.
247 Examples
248 --------
249 An example of the plot produced from this code is here:
251 .. image:: /_static/analysis_tools/wholeSkyPlotExample.png
253 For a detailed example of how to make a plot from the command line
254 please see the
255 :ref:`getting started guide<analysis-tools-getting-started>`.
256 """
257 skymap = kwargs["skymap"]
258 if plotInfo is None:
259 plotInfo = {}
261 if self.addThresholds:
262 metricThresholdFile = importResources.read_text(lsst.analysis.tools, "metricInformation.yaml")
263 metricDefs = yaml.safe_load(metricThresholdFile)
265 # Prevent Bands in the plot info showing a list of bands.
266 # If bands is a list, it implies that parameterizedBand=False,
267 # and that the metric is not band-specific.
268 if "bands" in plotInfo:
269 if isinstance(plotInfo["bands"], list):
270 plotInfo["bands"] = "N/A"
272 colorMap = self.colorMap
273 match self.colorMapType:
274 case "sequential":
275 if colorMap is None:
276 colorMap = stars_cmap()
277 outlierColor = "red"
278 case "divergent":
279 if colorMap is None:
280 colorMap = divergent_cmap()
281 outlierColor = "fuchsia"
283 # Setup figure.
284 if not self.publicationStyle:
285 figSize = (11, 3.5)
286 else:
287 figSize = (9, 3.5)
288 fig: AnnotatedFigure = make_figure(dpi=300, figsize=figSize)
289 set_rubin_plotstyle()
290 gs = gridspec.GridSpec(1, 20)
291 if not self.publicationStyle:
292 ax = fig.add_subplot(gs[3:15])
293 sp = skyproj.McBrydeSkyproj(ax=ax)
294 else:
295 ax = fig.add_subplot(111)
296 sp = skyproj.McBrydeSkyproj(ax=ax, extent=[-180, 180, -90, 30])
298 # This step is necessary to fix a race condition in skyproj
299 # It is fixed in a later version that is not yet in our env
300 proj = copy.copy(sp.ax.projection)
301 proj.set_plot_geodesics(False)
302 sp.draw_milky_way(label="Milky Way", zorder=101, alpha=0.3)
304 # Define color bar range.
305 colBarVals = data["z"]
306 if np.sum(np.isfinite(colBarVals)) > 0:
307 med = np.nanmedian(colBarVals)
308 else:
309 med = np.nan
310 sigmaMad = nanSigmaMad(colBarVals)
311 if self.colorBarMin is not None:
312 vmin = np.float64(self.colorBarMin)
313 else:
314 vmin = med - self.colorBarRange * sigmaMad
315 if self.colorBarMax is not None:
316 vmax = np.float64(self.colorBarMax)
317 else:
318 vmax = med + self.colorBarRange * sigmaMad
320 if self.fixAroundZero:
321 vmax = np.nanmax([np.fabs(vmin), np.fabs(vmax)])
322 vmin = -1 * vmax
324 dataName = self.zAxisLabel.format_map(kwargs)
325 if self.addThresholds and dataName in metricDefs:
326 if "lowThreshold" in metricDefs[dataName].keys():
327 lowThreshold = metricDefs[dataName]["lowThreshold"]
328 else:
329 lowThreshold = np.nan
330 if "highThreshold" in metricDefs[dataName].keys():
331 highThreshold = metricDefs[dataName]["highThreshold"]
332 else:
333 highThreshold = np.nan
334 outlierInds = np.where((colBarVals < lowThreshold) | (colBarVals > highThreshold))[0]
335 else:
336 # Note tracts with metrics outside (vmin, vmax) as outliers.
337 outlierInds = np.where((colBarVals < vmin) | (colBarVals > vmax))[0]
339 # Create patches using the corners of each tract.
340 tracts = []
341 ras = []
342 decs = []
343 mid_ras = []
344 mid_decs = []
345 allPatches = []
346 nanCount = 0
347 norm = Normalize(vmin=vmin, vmax=vmax, clip=True)
348 for i, tract in enumerate(data["tract"]):
349 corners = getTractCorners(skymap, tract)
350 tractRas = []
351 tractDecs = []
352 for ra, dec in corners:
353 tractRas.append(ra)
354 tractDecs.append(dec)
355 minRa = np.min(tractRas)
356 maxRa = np.max(tractRas)
357 minDec = np.min(tractDecs)
358 maxDec = np.max(tractDecs)
359 width = maxRa - minRa
360 height = maxDec - minDec
361 patch = Rectangle(
362 [minRa, minDec],
363 width,
364 height,
365 facecolor=colorMap(norm(data["z"][i])),
366 edgecolor="none",
367 zorder=20,
368 )
369 patch.set_transform(proj)
370 allPatches.append(patch)
371 sp.ax.add_patch(patch)
372 tracts.append(tract)
373 ras.append(corners[0][0])
374 decs.append(corners[0][1])
375 mid_ras.append((corners[0][0] + corners[1][0]) / 2)
376 mid_decs.append((corners[0][1] + corners[2][1]) / 2)
378 if self.showOutliers and not self.publicationStyle:
379 # Plot the outlier patches.
380 if i in outlierInds:
381 patch = Rectangle(
382 [minRa, minDec], width, height, facecolor="none", edgecolor=outlierColor, zorder=20
383 )
384 patch.set_transform(proj)
385 sp.ax.add_patch(patch)
387 if self.showNaNs:
388 # Plot tracts with NaN metric values
389 if not np.isfinite(colBarVals[i]):
390 nanCount += 1
391 patch = Rectangle(
392 [minRa, minDec],
393 width,
394 height,
395 facecolor="none",
396 edgecolor="grey",
397 zorder=20,
398 linewidth=0.5,
399 linestyle="dotted",
400 )
401 patch.set_transform(proj)
402 sp.ax.add_patch(patch)
404 # Add legend information.
405 if self.showOutliers and not self.publicationStyle:
406 sp.ax.plot([-1, -2], [-1, -2], color=outlierColor, label="Outlier")
407 if self.showNaNs:
408 sp.ax.plot([-1, -2], [-1, -2], color="grey", linestyle="dotted", label="NaN")
410 if self.labelTracts and not self.publicationStyle:
411 # Label the tracts
412 for i, tract in enumerate(tracts):
413 ax.text(
414 mid_ras[i],
415 mid_decs[i],
416 f"{tract}",
417 ha="center",
418 va="center",
419 fontsize=2,
420 alpha=0.7,
421 zorder=100,
422 )
424 if not self.publicationStyle:
425 ax1 = fig.add_axes([0.73, 0.25, 0.20, 0.47])
427 if np.sum(np.isfinite(data["z"])) > 0:
428 ax1.hist(data["z"], bins=len(data["z"] / 10), color=stars_color(), histtype="step")
429 else:
430 ax1.text(0.5, 0.5, "Data all NaN/Inf")
431 ax1.set_xlabel("Metric Values")
432 ax1.set_ylabel("Number")
433 ax1.yaxis.set_label_position("right")
434 ax1.yaxis.tick_right()
436 if self.addThresholds and dataName in metricDefs and not self.publicationStyle:
437 # Check the thresholds are finite and set them to
438 # the min/max of the data if they aren't to calculate
439 # the x range of the plot
440 if np.isfinite(lowThreshold):
441 ax1.axvline(lowThreshold, color=accent_color())
442 else:
443 lowThreshold = np.nanmin(colBarVals)
444 if np.isfinite(highThreshold):
445 ax1.axvline(highThreshold, color=accent_color())
446 else:
447 highThreshold = np.nanmax(colBarVals)
449 widthThreshold = highThreshold - lowThreshold
450 upperLim = highThreshold + 0.5 * widthThreshold
451 lowerLim = lowThreshold - 0.5 * widthThreshold
452 ax1.set_xlim(lowerLim, upperLim)
453 numOutside = np.sum((data["z"] > upperLim) | (data["z"] < lowerLim))
454 ax1.set_title("Outside plot limits: " + str(numOutside))
456 else:
457 if vmin != vmax and np.isfinite(vmin) and np.isfinite(vmax):
458 ax1.set_xlim(vmin, vmax)
460 if self.autoAxesLimits:
461 xlim, ylim = self._getAxesLimits(ras, decs)
462 else:
463 xlim, ylim = self.xLimits, self.yLimits
465 sp.set_extent([xlim[0], xlim[1], ylim[0], ylim[1]])
467 if self.showOutliers and not self.publicationStyle:
468 # Add text boxes to show the number of tracts, number of NaNs,
469 # median, sigma MAD, and the five largest outlier values.
470 outlierText = self._getMaxOutlierVals(self.colorBarRange, tracts, colBarVals, outlierInds)
471 if not self.publicationStyle:
472 # Make vertical text spacing readable for different figure sizes.
473 multiplier = 3.5 / fig.get_size_inches()[1]
474 verticalSpacing = 0.028 * multiplier
475 fig.text(
476 0.01,
477 0.01 + 3 * verticalSpacing,
478 f"Num tracts: {len(tracts)}",
479 transform=fig.transFigure,
480 fontsize=8,
481 alpha=0.7,
482 )
484 if self.showNaNs and not self.publicationStyle:
485 fig.text(
486 0.01,
487 0.01 + 2 * verticalSpacing,
488 f"Num nans: {nanCount}",
489 transform=fig.transFigure,
490 fontsize=8,
491 alpha=0.7,
492 )
493 if not self.publicationStyle:
494 fig.text(
495 0.01,
496 0.01 + verticalSpacing,
497 f"Median: {med:.3f}; " + r"$\sigma_{MAD}$" + f": {sigmaMad:.3f}",
498 transform=fig.transFigure,
499 fontsize=8,
500 alpha=0.7,
501 )
502 if self.showOutliers and not self.publicationStyle:
503 fig.text(0.01, 0.01, outlierText, transform=fig.transFigure, fontsize=8, alpha=0.7)
505 # Make the color bar with a metric label.
506 cax = fig.add_axes([0.04, 0.15, 0.03, 0.7])
507 cr = colorizer.Colorizer(norm=norm, cmap=colorMap)
508 fig.colorbar(
509 colorizer.ColorizingArtist(cr),
510 cax=cax,
511 shrink=0.7,
512 extend="both",
513 orientation="vertical",
514 )
515 cbarText = "Metric Values"
516 cax.tick_params(labelsize=8)
518 text = cax.text(
519 0.5,
520 0.5,
521 cbarText,
522 transform=cax.transAxes,
523 ha="center",
524 va="center",
525 fontsize=10,
526 zorder=100,
527 rotation="vertical",
528 )
529 text.set_path_effects([pathEffects.Stroke(linewidth=3, foreground="w"), pathEffects.Normal()])
531 # Finalize plot appearance.
532 if not self.publicationStyle:
533 addPlotInfo(fig, plotInfo)
534 fig.subplots_adjust(left=0.02, right=0.92, top=0.8, bottom=0.17, wspace=0.1)
535 else:
536 fig.subplots_adjust(left=0.2, right=0.85, top=0.8, bottom=0.17, wspace=0.1)
537 titleText = self.zAxisLabel.format_map(kwargs)
538 if "zUnit" in data and data["zUnit"] != "":
539 titleText += f" ({data['zUnit']})"
540 fig.suptitle("Metric: " + titleText, fontsize=12)
542 if self.publicationStyle:
543 sp.ax.legend(loc="upper left", bbox_to_anchor=(0.9, 0.1))
544 labelFontsize = 15
545 else:
546 sp.ax.legend(loc="upper right", bbox_to_anchor=(1.47, 1.56), ncols=2)
547 labelFontsize = 15
548 sp.ax.set_xlabel(self.xAxisLabel, fontsize=labelFontsize)
549 sp.ax.set_ylabel(self.yAxisLabel, fontsize=labelFontsize, labelpad=10)
550 sp.ax.tick_params(labelsize=8)
551 sp.ax.tick_params(axis="y", labelsize=8)
553 # This saves metadata in the PNG that allows the plot-navigator
554 # to provide tract numbers and metric values on mouseover.
555 #
556 # PNG metadata is a set of string keys and string values.
557 # The WholeSkyPlot stores two keys:
558 # - label: the string describing the regions ('tract')
559 # - boxes, JSON string of a list of per-region dictionaries,
560 # where each dictionary has fields:
561 # - min_x, max_x, min_y, max_y for the pixel coordinates of
562 # the four corners of the region
563 # - id: the identifier of the region (e.g. tract number)
564 # - value: the region's metric, as a string.
565 #
566 def make_patch_md(patch, id_field, value, ax):
567 trans = patch.get_transform()
568 path = trans.transform_path(patch.get_path())
569 scale = fig.canvas.device_pixel_ratio
570 x_path = [int(x) for x in (path.vertices[:, 0] / scale).tolist()]
571 y_path = [int(y) for y in (path.vertices[:, 1] / scale).tolist()]
572 return {
573 "min_x": min(x_path),
574 "max_x": max(x_path),
575 "min_y": min(y_path),
576 "max_y": max(y_path),
577 "id": f"{id_field}",
578 "value": f"{value:.3}",
579 }
581 # After ax.set_aspect(), the figure needs to be drawn for the axes
582 # transformations to be updated to the right values.
583 fig.canvas.draw_idle()
585 patch_coordinate_entries = [
586 make_patch_md(patch, tract, value, ax)
587 for (patch, tract, value) in zip(allPatches, tracts, colBarVals)
588 ]
590 fig.metadata = {"label": "Tract", "boxes": json.dumps(patch_coordinate_entries)}
592 return fig