Coverage for python / lsst / analysis / tools / actions / plot / skyPlot.py: 11%
197 statements
« prev ^ index » next coverage.py v7.13.2, created at 2026-02-19 01:44 -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__ = ("SkyPlot",)
26from collections.abc import Mapping
28import matplotlib.patheffects as pathEffects
29import numpy as np
30from matplotlib.figure import Figure
31from matplotlib.patches import Rectangle
33from lsst.pex.config import Field, ListField
34from lsst.pex.config.configurableActions import ConfigurableActionField
35from lsst.utils.plotting import (
36 divergent_cmap,
37 galaxies_cmap,
38 galaxies_color,
39 make_figure,
40 set_rubin_plotstyle,
41 stars_cmap,
42 stars_color,
43)
45from ...interfaces import KeyedData, KeyedDataSchema, PlotAction, Scalar, Vector, VectorAction
46from ...math import nanMedian, nanSigmaMad
47from .calculateRange import Med2Mad
48from .plotUtils import addPlotInfo, generateSummaryStats, plotProjectionWithBinning, sortAllArrays
51class SkyPlot(PlotAction):
52 """Plots the on sky distribution of a parameter.
54 Plots the values of the parameter given for the z axis
55 according to the positions given for x and y. Optimised
56 for use with RA and Dec. Also calculates some basic
57 statistics and includes those on the plot.
59 The plotting of patch outlines requires patch information
60 to be included as an additional parameter.
61 """
63 xAxisLabel = Field[str](doc="Label to use for the x axis.", optional=False)
64 yAxisLabel = Field[str](doc="Label to use for the y axis.", optional=False)
65 zAxisLabel = Field[str](doc="Label to use for the z axis.", optional=False)
67 plotOutlines = Field[bool](
68 doc="Plot the outlines of the ccds/patches?",
69 default=True,
70 )
72 plotTypes = ListField[str](
73 doc="Selection of types of objects to plot. Can take any combination of"
74 " stars, galaxies, unknown, mag, any.",
75 optional=False,
76 # itemCheck=_validatePlotTypes,
77 )
79 plotName = Field[str](doc="The name for the plot.", optional=False)
81 fixAroundZero = Field[bool](
82 doc="Fix the colorbar to be symmetric around zero.",
83 default=False,
84 )
86 colorbarRange = ConfigurableActionField[VectorAction](
87 doc="Action to calculate the min and max of the colorbar range.",
88 default=Med2Mad,
89 )
91 showExtremeOutliers = Field[bool](
92 doc="Show the x-y positions of extreme outlier values as overlaid scatter points.",
93 default=True,
94 )
96 publicationStyle = Field[bool](
97 doc="Make a simplified plot for publication use.",
98 default=False,
99 )
101 divergent = Field[bool](
102 doc="Use a divergent colormap?",
103 default=False,
104 )
106 def getInputSchema(self, **kwargs) -> KeyedDataSchema:
107 base = []
108 if "stars" in self.plotTypes: # type: ignore
109 base.append(("xStars", Vector))
110 base.append(("yStars", Vector))
111 base.append(("zStars", Vector))
112 base.append(("starStatMask", Vector))
113 if "galaxies" in self.plotTypes: # type: ignore
114 base.append(("xGalaxies", Vector))
115 base.append(("yGalaxies", Vector))
116 base.append(("zGalaxies", Vector))
117 base.append(("galaxyStatMask", Vector))
118 if "unknown" in self.plotTypes: # type: ignore
119 base.append(("xUnknowns", Vector))
120 base.append(("yUnknowns", Vector))
121 base.append(("zUnknowns", Vector))
122 base.append(("unknownStatMask", Vector))
123 if "any" in self.plotTypes: # type: ignore
124 base.append(("x", Vector))
125 base.append(("y", Vector))
126 base.append(("z", Vector))
127 base.append(("statMask", Vector))
129 return base
131 def __call__(self, data: KeyedData, **kwargs) -> Mapping[str, Figure] | Figure:
132 self._validateInput(data, **kwargs)
133 return self.makePlot(data, **kwargs)
134 # table is a dict that needs: x, y, run, skymap, filter, tract,
136 def _validateInput(self, data: KeyedData, **kwargs) -> None:
137 """NOTE currently can only check that something is not a Scalar, not
138 check that the data is consistent with Vector
139 """
140 needed = self.getInputSchema(**kwargs)
141 if remainder := {key.format(**kwargs) for key, _ in needed} - {
142 key.format(**kwargs) for key in data.keys()
143 }:
144 raise ValueError(f"Task needs keys {remainder} but they were not found in input")
145 for name, typ in needed:
146 isScalar = issubclass((colType := type(data[name.format(**kwargs)])), Scalar)
147 if isScalar and typ != Scalar:
148 raise ValueError(f"Data keyed by {name} has type {colType} but action requires type {typ}")
150 def statsAndText(self, arr, mask=None):
151 """Calculate some stats from an array and return them
152 and some text.
153 """
154 numPoints = len(arr)
155 if mask is not None:
156 arr = arr[mask]
157 med = nanMedian(arr)
158 sigMad = nanSigmaMad(arr)
160 statsText = (
161 f"Median: {med:0.2f}\n"
162 + r"$\sigma_{MAD}$: "
163 + f"{sigMad:0.2f}\n"
164 + r"n$_{points}$: "
165 + f"{numPoints}"
166 )
168 return med, sigMad, statsText
170 def makePlot(
171 self,
172 data: KeyedData,
173 plotInfo: Mapping[str, str] | None = None,
174 sumStats: Mapping | None = None,
175 **kwargs,
176 ) -> Figure:
177 """Make a skyPlot of the given data.
179 Parameters
180 ----------
181 data : `KeyedData`
182 The catalog to plot the points from.
183 plotInfo : `dict`
184 A dictionary of information about the data being plotted with keys:
186 ``"run"``
187 The output run for the plots (`str`).
188 ``"skymap"``
189 The type of skymap used for the data (`str`).
190 ``"filter"``
191 The filter used for this data (`str`).
192 ``"tract"``
193 The tract that the data comes from (`str`).
195 sumStats : `dict`
196 A dictionary where the patchIds are the keys which store the R.A.
197 and dec of the corners of the patch.
199 Returns
200 -------
201 `pipeBase.Struct` containing:
202 skyPlot : `matplotlib.figure.Figure`
203 The resulting figure.
205 Notes
206 -----
207 Expects the data to contain slightly different things
208 depending on the types specified in plotTypes. This
209 is handled automatically if you go through the pipetask
210 framework but if you call this method separately then you
211 need to make sure that data contains what the code is expecting.
213 If stars is in the plot types given then it is expected that
214 data contains: xStars, yStars, zStars and starStatMask.
216 If galaxies is present: xGalaxies, yGalaxies, zGalaxies and
217 galaxyStatsMask.
219 If unknown is present: xUnknowns, yUnknowns, zUnknowns and
220 unknownStatMask.
222 If any is specified: x, y, z, statMask.
224 These options are not exclusive and multiple can be specified
225 and thus need to be present in data.
227 Examples
228 --------
229 An example of the plot produced from this code is here:
231 .. image:: /_static/analysis_tools/skyPlotExample.png
233 For a detailed example of how to make a plot from the command line
234 please see the
235 :ref:`getting started guide<analysis-tools-getting-started>`.
236 """
238 set_rubin_plotstyle()
239 fig = make_figure()
240 ax = fig.add_subplot(111)
242 if sumStats is None:
243 if self.plotOutlines and "patch" in data.keys():
244 sumStats = generateSummaryStats(data, kwargs["skymap"], plotInfo)
245 else:
246 sumStats = {}
248 if plotInfo is None:
249 plotInfo = {}
251 # Make divergent colormaps for stars, galaxes and all the points
252 starsCmap = stars_cmap()
253 galsCmap = galaxies_cmap()
255 xCol = self.xAxisLabel
256 yCol = self.yAxisLabel
257 zCol = self.zAxisLabel # noqa: F841
259 toPlotList = []
260 # For galaxies
261 if "galaxies" in self.plotTypes:
262 sortedArrs = sortAllArrays(
263 [data["zGalaxies"], data["xGalaxies"], data["yGalaxies"], data["galaxyStatMask"]]
264 )
265 [colorValsGalaxies, xsGalaxies, ysGalaxies, statGalaxies] = sortedArrs
266 statGalMed, statGalMad, galStatsText = self.statsAndText(colorValsGalaxies, mask=statGalaxies)
267 # Add statistics
268 bbox = dict(facecolor=galaxies_color(), alpha=0.5, edgecolor="none")
269 # Check if plotting stars and galaxies, if so move the
270 # text box so that both can be seen. Needs to be
271 # > 2 becuase not being plotted points are assigned 0
272 if len(self.plotTypes) > 2:
273 boxLoc = 0.63
274 else:
275 boxLoc = 0.8
276 ax.text(boxLoc, 0.91, galStatsText, transform=fig.transFigure, fontsize=8, bbox=bbox)
277 if self.divergent:
278 galsCmap = divergent_cmap()
279 toPlotList.append((xsGalaxies, ysGalaxies, colorValsGalaxies, galsCmap, "Galaxies"))
281 # For stars
282 if "stars" in self.plotTypes:
283 sortedArrs = sortAllArrays([data["zStars"], data["xStars"], data["yStars"], data["starStatMask"]])
284 [colorValsStars, xsStars, ysStars, statStars] = sortedArrs
285 statStarMed, statStarMad, starStatsText = self.statsAndText(colorValsStars, mask=statStars)
286 if not self.publicationStyle:
287 # Add statistics
288 bbox = dict(facecolor=stars_color(), alpha=0.5, edgecolor="none")
289 ax.text(0.8, 0.91, starStatsText, transform=fig.transFigure, fontsize=8, bbox=bbox)
290 if self.divergent:
291 starsCmap = divergent_cmap()
292 toPlotList.append((xsStars, ysStars, colorValsStars, starsCmap, "Stars"))
294 # For unknowns
295 if "unknown" in self.plotTypes:
296 sortedArrs = sortAllArrays(
297 [data["zUnknowns"], data["xUnknowns"], data["yUnknowns"], data["unknownStatMask"]]
298 )
299 [colorValsUnknowns, xsUnknowns, ysUnknowns, statUnknowns] = sortedArrs
300 statUnknownMed, statUnknownMad, unknownStatsText = self.statsAndText(
301 colorValsUnknowns, mask=statUnknowns
302 )
303 if not self.publicationStyle:
304 bbox = dict(facecolor="green", alpha=0.2, edgecolor="none")
305 ax.text(0.8, 0.91, unknownStatsText, transform=fig.transFigure, fontsize=8, bbox=bbox)
306 toPlotList.append((xsUnknowns, ysUnknowns, colorValsUnknowns, "viridis", "Unknown"))
308 if "any" in self.plotTypes:
309 sortedArrs = sortAllArrays([data["z"], data["x"], data["y"], data["statMask"]])
310 [colorValsAny, xs, ys, statAny] = sortedArrs
311 statAnyMed, statAnyMad, anyStatsText = self.statsAndText(colorValsAny, mask=statAny)
312 if not self.publicationStyle:
313 bbox = dict(facecolor="#bab0ac", alpha=0.2, edgecolor="none")
314 ax.text(0.8, 0.91, anyStatsText, transform=fig.transFigure, fontsize=8, bbox=bbox)
315 toPlotList.append((xs, ys, colorValsAny, "viridis", "All"))
317 # Corner plot of patches showing summary stat in each
318 if self.plotOutlines:
319 patches = []
320 for dataId in sumStats.keys():
321 (corners, _) = sumStats[dataId]
322 ra = corners[0][0].asDegrees()
323 dec = corners[0][1].asDegrees()
324 xy = (ra, dec)
325 width = corners[2][0].asDegrees() - ra
326 height = corners[2][1].asDegrees() - dec
327 patches.append(Rectangle(xy, width, height, alpha=0.3))
328 ras = [ra.asDegrees() for (ra, dec) in corners]
329 decs = [dec.asDegrees() for (ra, dec) in corners]
330 ax.plot(ras + [ras[0]], decs + [decs[0]], "k", lw=0.5)
331 cenX = ra + width / 2
332 cenY = dec + height / 2
333 if dataId == "tract":
334 minRa = np.min(ras)
335 minDec = np.min(decs)
336 maxRa = np.max(ras)
337 maxDec = np.max(decs)
338 if dataId != "tract":
339 ax.annotate(
340 dataId,
341 (cenX, cenY),
342 color="k",
343 fontsize=5,
344 ha="center",
345 va="center",
346 path_effects=[pathEffects.withStroke(linewidth=2, foreground="w")],
347 )
349 for i, (xs, ys, colorVals, cmap, label) in enumerate(toPlotList):
350 finite = np.isfinite(xs) & np.isfinite(ys)
351 xs = xs[finite]
352 ys = ys[finite]
353 colorVals = colorVals[finite]
354 n_xs = len(xs)
355 if n_xs == 0:
356 continue
358 # colorVal column is unusable so zero it out
359 # This should be obvious on the plot
360 if not any(np.isfinite(colorVals)):
361 colorVals[:] = 0
362 minColorVal, maxColorVal = self.colorbarRange(colorVals)
364 if not self.plotOutlines or "tract" not in sumStats.keys():
365 minRa = np.min(xs)
366 maxRa = np.max(xs)
367 minDec = np.min(ys)
368 maxDec = np.max(ys)
369 # Avoid identical end points which causes problems in binning
370 if minRa == maxRa:
371 maxRa += 1e-5 # There is no reason to pick this number in particular
372 if minDec == maxDec:
373 maxDec += 1e-5 # There is no reason to pick this number in particular
374 if n_xs < 5:
375 continue
377 if self.publicationStyle:
378 showExtremeOutliers = False
379 else:
380 showExtremeOutliers = self.showExtremeOutliers
381 plotOut = plotProjectionWithBinning(
382 ax,
383 xs,
384 ys,
385 colorVals,
386 cmap,
387 minRa,
388 maxRa,
389 minDec,
390 maxDec,
391 vmin=minColorVal,
392 vmax=maxColorVal,
393 fixAroundZero=self.fixAroundZero,
394 isSorted=True,
395 showExtremeOutliers=showExtremeOutliers,
396 )
397 ax.set_aspect("equal")
398 if not self.publicationStyle:
399 cax = fig.add_axes([0.87 + i * 0.04, 0.11, 0.04, 0.77])
400 fig.colorbar(plotOut, cax=cax, extend="both")
401 else:
402 fig.subplots_adjust(wspace=0.0, hspace=0.0, right=0.95, bottom=0.15)
403 axBbox = ax.get_position()
404 cax = fig.add_axes([axBbox.x1, axBbox.y0, 0.04, axBbox.y1 - axBbox.y0])
405 fig.colorbar(plotOut, cax=cax)
407 colorBarLabel = f"{self.zAxisLabel}: {label}"
408 text = cax.text(
409 0.5,
410 0.5,
411 colorBarLabel,
412 color="k",
413 rotation="vertical",
414 transform=cax.transAxes,
415 ha="center",
416 va="center",
417 fontsize=10,
418 )
419 text.set_path_effects([pathEffects.Stroke(linewidth=3, foreground="w"), pathEffects.Normal()])
421 if i == 0 and len(toPlotList) > 1:
422 cax.yaxis.set_ticks_position("left")
424 ax.set_xlabel(xCol)
425 ax.set_ylabel(yCol)
427 fig.canvas.draw()
429 # Find some useful axis limits
430 lenXs = [len(xs) for (xs, _, _, _, _) in toPlotList]
431 if lenXs != [] and np.max(lenXs) > 1000:
432 padRa = (maxRa - minRa) / 10
433 padDec = (maxDec - minDec) / 10
434 ax.set_xlim(maxRa + padRa, minRa - padRa)
435 ax.set_ylim(minDec - padDec, maxDec + padDec)
436 else:
437 ax.invert_xaxis()
439 # Add useful information to the plot
440 if not self.publicationStyle:
441 fig.subplots_adjust(wspace=0.0, hspace=0.0)
442 fig = addPlotInfo(fig, plotInfo)
444 return fig