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