Coverage for python/lsst/analysis/ap/spatiallySampledMetricsQA.py: 10%
138 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-09-24 02:27 -0700
1# This file is part of analysis_ap.
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/>.
22"""Notebook-friendly QA for ``SpatiallySampledMetricsTask`` output.
24The single public entry point `subtraction_quality_report` consumes a
25``*_spatiallySampledMetrics`` table for one detector and produces:
27- A printed top-line summary of the three diagnostic scalars
28 (`diffim_chi2PerPix`, `psfMatchingKernel_residualNorm`,
29 `dipole_density`) at percentiles that make localized failures visible.
30- A three-panel diagnostic figure with the three scalars rendered as
31 linear-interpolated heatmaps, with the kernel centroid offset quiver
32 (colored cyclically by angle, magnitude shown by a per-panel
33 quiverkey reference arrow) overlaid on each panel.
35See the module-level constants for the metric reference values and the
36default mask-fraction columns used to filter samples sitting on bad
37detector regions before computing statistics.
38"""
40from __future__ import annotations
42__all__ = ["subtraction_quality_report"]
44import numpy as np
45import pandas as pd
47from lsst.analysis.ap.skymapOverlay import make_affine_sky_to_xy, draw_skymap_outlines_mpl
49# Mask-fraction columns whose sum indicates a sample sits on an unusable
50# region of the detector. The headline scalars are computed only on samples
51# below ``bad_mask_threshold`` so the distribution tails reflect subtraction
52# quality rather than edge / saturated pixels.
53DEFAULT_BAD_MASK_COLUMNS = (
54 "bad_mask_fraction",
55 "sat_mask_fraction",
56 "edge_mask_fraction",
57 "no_data_mask_fraction",
58)
60# Default extra padding (in detector pixels) added to each panel beyond
61# the autoscaled data limits, so labels anchored at the edge of the data have
62# room to render without being clipped by the panel boundary. NOTE: the
63# displayed figure is at significantly lower resolution than the original
64# image.
65_DEFAULT_PANEL_PADDING_PIX = 150
67# (column, reference value, display label). Reference is the value the
68# metric takes on a perfectly subtracted, well-decorrelated diffim.
69_HEADLINE_METRICS = (
70 ("diffim_chi2PerPix", 1.0, "Diffim chi^2/pix"), # noqa: E241
71 ("psfMatchingKernel_residualNorm", 0.0, "PSF match residual"), # noqa: E241
72 ("dipole_density", 0.0, "Dipoles / deg^2"), # noqa: E241
73)
76def _coerce_to_frame(metrics):
77 """Accept an astropy Table, DataFrame, or any DataFrame-castable input."""
78 if isinstance(metrics, pd.DataFrame):
79 return metrics
80 if hasattr(metrics, "to_pandas"):
81 return metrics.to_pandas()
82 return pd.DataFrame(metrics)
85def _filter_clean(df, bad_mask_threshold, bad_mask_columns):
86 """Drop samples whose summed bad-mask fractions exceed the threshold."""
87 cols = [c for c in bad_mask_columns if c in df.columns]
88 if not cols:
89 return df.copy()
90 return df[df[cols].sum(axis=1) < bad_mask_threshold].copy()
93def _percentile_row(values):
94 """Return formatted (median, p84, p95, p99) strings for a 1-D array."""
95 finite = values[np.isfinite(values)]
96 if finite.size == 0:
97 return ["nan"]*4
98 q = np.percentile(finite, [50, 84, 95, 99])
99 return [f"{x:.3f}" for x in q]
102def _print_summary(clean, n_total, threshold):
103 """Print a fixed-width summary of the three headline metrics."""
104 headers = ["metric", "ref", "median", "p84", "p95", "p99"]
105 rows = []
106 for col, ref, _label in _HEADLINE_METRICS:
107 if col not in clean.columns:
108 rows.append([col, f"{ref:.2f}", "n/a", "n/a", "n/a", "n/a"])
109 continue
110 rows.append([col, f"{ref:.2f}", *_percentile_row(clean[col].to_numpy())])
111 widths = [max(len(str(r[i])) for r in [headers, *rows]) for i in range(len(headers))]
113 def _fmt(row):
114 return " ".join(str(c).ljust(w) for c, w in zip(row, widths))
116 print(f"SpatiallySampledMetrics: {len(clean)}/{n_total} samples retained "
117 f"(bad-mask fraction sum < {threshold:g})")
118 print()
119 print(_fmt(headers))
120 print(" ".join("-"*w for w in widths))
121 for row in rows:
122 print(_fmt(row))
125def _metric_panel(ax, fig, clean, col, vmin, vmax, label, cmap, vcenter=None,
126 panel_padding_pix=_DEFAULT_PANEL_PADDING_PIX):
127 """Render one metric panel: linear-interpolated heatmap with sample
128 markers.
130 Points whose (x, y, value) tuple contains any NaN are dropped before
131 interpolation. Grid cells outside the convex hull of the surviving
132 samples are left as NaN, which the colormap renders as transparent.
134 If ``vcenter`` is provided and strictly between ``vmin`` and ``vmax``,
135 the panel uses a ``TwoSlopeNorm`` so the colormap midpoint always
136 maps to ``vcenter`` regardless of whether the (vmin, vmax) range is
137 symmetric about it.
138 """
139 if col not in clean.columns or not clean[col].notna().any():
140 ax.text(0.5, 0.5, f"{col}\n(not present)",
141 ha="center", va="center", transform=ax.transAxes)
142 ax.set_title(col)
143 return
145 x = clean["x"].to_numpy()
146 y = clean["y"].to_numpy()
147 z = clean[col].to_numpy()
148 valid = np.isfinite(x) & np.isfinite(y) & np.isfinite(z)
149 if not valid.any():
150 ax.text(0.5, 0.5, f"{col}\n(no valid samples)",
151 ha="center", va="center", transform=ax.transAxes)
152 ax.set_title(col)
153 return
154 xv, yv, zv = x[valid], y[valid], z[valid]
156 # vmin and vmax autoscale independently: an explicit ``vmin=0.0`` from
157 # the caller is never overridden by the vmax-None branch.
158 if vmin is None:
159 vmin = 0.0
160 if vmax is None:
161 vmax = float(np.nanpercentile(zv, 99))
162 # Set a reasonable colorbar scale even if the data is completely constant.
163 if not (vmax > vmin):
164 vmax = vmin + 1.0
166 from scipy.interpolate import griddata
167 from matplotlib.colors import Normalize, TwoSlopeNorm
168 grid_size = 200
169 xi = np.linspace(xv.min(), xv.max(), grid_size)
170 yi = np.linspace(yv.min(), yv.max(), grid_size)
171 XI, YI = np.meshgrid(xi, yi)
172 ZI = griddata((xv, yv), zv, (XI, YI), method="linear")
174 # Center the colormap on ``vcenter`` when supplied and well-posed.
175 # TwoSlopeNorm requires vmin < vcenter < vmax strictly; if the
176 # autoscaled range collapses around vcenter we fall back to a plain
177 # clipping Normalize rather than raise.
178 if vcenter is not None and vmin < vcenter < vmax:
179 norm = TwoSlopeNorm(vcenter=vcenter, vmin=vmin, vmax=vmax)
180 else:
181 norm = Normalize(vmin=vmin, vmax=vmax, clip=True)
182 im = ax.imshow(ZI, origin="lower",
183 extent=(xv.min(), xv.max(), yv.min(), yv.max()),
184 norm=norm, cmap=cmap, aspect="equal")
185 fig.colorbar(im, ax=ax, label=label)
186 # Overlay the kernel centroid quiver so direction information is
187 # available next to every metric heatmap, with a quiverkey reference
188 # arrow. Sample markers go on top with higher "zorder" so the positions
189 # stay visible through arrows.
190 quiver_info = _overlay_kernel_quiver(ax, clean)
191 if quiver_info is not None:
192 q, ref_length = quiver_info
193 ax.quiverkey(q, 0.85, 1.05, ref_length, f"{ref_length:.2f}″",
194 labelpos="E", coordinates="axes")
195 ax.scatter(xv, yv, s=8, facecolors="white", edgecolors="black",
196 linewidths=0.3, zorder=5)
197 ax.set_title(col)
198 ax.set_xlabel("x [pix]")
199 ax.set_ylabel("y [pix]")
200 ax.set_aspect("equal")
201 # Pad the data limits relative to whatever the artists left them at,
202 # so labels (and quiver heads) sitting at the very edge of the data
203 # are not clipped by the panel boundary.
204 xmin, xmax = ax.get_xlim()
205 ymin, ymax = ax.get_ylim()
206 ax.set_xlim(xmin - panel_padding_pix, xmax + panel_padding_pix)
207 ax.set_ylim(ymin - panel_padding_pix, ymax + panel_padding_pix)
210def _overlay_kernel_quiver(ax, clean):
211 """Draw the HSV-colored kernel centroid quiver onto ``ax``.
213 Returns
214 -------
215 info : tuple or None
216 ``(quiver_artist, ref_length_arcsec)`` on success, or None if the
217 required columns are missing or no samples are valid. The caller
218 decides whether to render a colorbar / quiverkey for it.
219 """
220 needed = ("psfMatchingKernel_length", "psfMatchingKernel_direction", "x", "y")
221 if not all(c in clean.columns for c in needed):
222 return None
223 length = clean["psfMatchingKernel_length"].to_numpy()
224 direction = clean["psfMatchingKernel_direction"].to_numpy()
225 u = length*np.cos(direction)
226 v = length*np.sin(direction)
227 ok = np.isfinite(u) & np.isfinite(v) & (length > 0)
228 if not ok.any():
229 return None
231 ref_length = float(np.nanpercentile(length[ok], 95))
232 direction_deg = np.degrees(direction[ok]) % 360.0
233 x = clean["x"].to_numpy()[ok]
234 y = clean["y"].to_numpy()[ok]
235 q = ax.quiver(x, y, u[ok], v[ok], direction_deg,
236 cmap="hsv", clim=(0, 360),
237 angles="xy", scale_units="xy",
238 scale=ref_length/200, width=0.004, pivot="mid",
239 zorder=4)
240 return q, ref_length
243def _draw_skymap_outlines(ax, skymap, clean, label_fontsize=7):
244 """Overlay patch boundaries (with tract,patch labels) on a panel.
246 The metrics table carries no WCS, so the sky↔detector pixel mapping is
247 derived from the sample positions' ``(x, y)`` and
248 ``(coord_ra, coord_dec)`` columns via a least-squares affine fit (see
249 `~lsst.analysis.ap.skymapOverlay.make_affine_sky_to_xy`). For a single
250 detector this is typically accurate to well under a pixel — enough for
251 visualization but not for science.
253 Silently no-ops when sky coordinates aren't available or fewer than
254 three samples are valid.
255 """
256 if "coord_ra" not in clean.columns or "coord_dec" not in clean.columns:
257 return
259 ra = clean["coord_ra"].to_numpy()
260 dec = clean["coord_dec"].to_numpy()
261 xs = clean["x"].to_numpy()
262 ys = clean["y"].to_numpy()
263 valid = np.isfinite(ra) & np.isfinite(dec) & np.isfinite(xs) & np.isfinite(ys)
264 if valid.sum() < 3:
265 return
267 import lsst.geom as geom
269 sky_to_xy = make_affine_sky_to_xy(ra[valid], dec[valid], xs[valid], ys[valid])
271 # Build a coord list spanning the sample footprint so findTractPatchList
272 # returns every tract / patch that touches it.
273 corner_pairs = [
274 (ra[valid].min(), dec[valid].min()),
275 (ra[valid].max(), dec[valid].min()),
276 (ra[valid].max(), dec[valid].max()),
277 (ra[valid].min(), dec[valid].max()),
278 ]
279 sky_corners = [geom.SpherePoint(r, d, geom.radians) for r, d in corner_pairs]
281 draw_skymap_outlines_mpl(ax, skymap, sky_to_xy, sky_corners,
282 label_fontsize=label_fontsize)
285def _make_figure(clean, chi2_scale=1.0, skymap=None, label_fontsize=7,
286 panel_padding_pix=_DEFAULT_PANEL_PADDING_PIX):
287 """Build a three-panel diagnostic figure.
289 Parameters
290 ----------
291 clean : `pandas.DataFrame`
292 The filtered metrics table.
293 chi2_scale : `float`, optional
294 Half-width of the ``diffim_chi2PerPix`` colormap range, measured
295 multiplicatively about the nominal value of 1.0. ``vmax`` is set
296 to ``1 + chi2_scale`` and ``vmin`` to its reciprocal so the
297 colorbar covers the same factor above and below nominal in log space.
298 skymap : `lsst.skymap.BaseSkyMap`, optional
299 If supplied, overlay patch outlines (with tract,patch labels) on
300 every panel.
301 label_fontsize : `int` or `float`, optional
302 Font size for the per-patch ``tract,patch`` labels.
303 """
304 import matplotlib.pyplot as plt
305 from matplotlib.colors import LinearSegmentedColormap
307 # White-to-black sequential colormap for dipole density: pure white
308 # at the floor, pure black at vmax. NaN cells (outside the
309 # interpolation convex hull) are mapped to ``lightblue`` instead of
310 # any grey because every grey lives somewhere inside the white→black
311 # gradient and would otherwise read as a real mid-range value.
312 white_to_black = LinearSegmentedColormap.from_list(
313 "white_to_black", [(1.0, 1.0, 1.0), (0.0, 0.0, 0.0)]).copy()
314 white_to_black.set_bad("lightblue")
316 chi2_vmax = 1.0 + chi2_scale
317 chi2_vmin = 1.0/chi2_vmax
319 fig, axes = plt.subplots(1, 3, figsize=(15, 4.5), constrained_layout=True)
320 # (col, vmin, vmax, label, cmap, vcenter). ``vcenter`` is non-None only
321 # for panels whose colormap should be anchored at a known reference
322 # value -- e.g. chi^2/pix is centered on its nominal value of 1.0.
323 panels = (
324 ("diffim_chi2PerPix", chi2_vmin, chi2_vmax, "Diffim chi^2/pix", "RdBu_r", 1.0), # noqa: E241,E501
325 ("psfMatchingKernel_residualNorm", 0.0, 0.2, "PSF match residual", "viridis", None), # noqa: E241,E501
326 ("dipole_density", 0.0, None, "Dipoles / deg^2", white_to_black, None), # noqa: E241,E501
327 )
328 for ax, (col, vmin, vmax, label, cmap, vcenter) in zip(axes, panels):
329 _metric_panel(ax, fig, clean, col, vmin, vmax, label, cmap, vcenter=vcenter,
330 panel_padding_pix=panel_padding_pix)
331 if skymap is not None:
332 _draw_skymap_outlines(ax, skymap, clean, label_fontsize=label_fontsize)
333 return fig
336def subtraction_quality_report(metrics,
337 bad_mask_threshold=0.2,
338 bad_mask_columns=DEFAULT_BAD_MASK_COLUMNS,
339 chi2_scale=1.0,
340 skymap=None,
341 label_fontsize=7,
342 panel_padding_pix=_DEFAULT_PANEL_PADDING_PIX):
343 """Print a headline metric summary and build the diagnostic plot.
345 Parameters
346 ----------
347 metrics : `astropy.table.Table` or `pandas.DataFrame`
348 The metrics table produced by ``SpatiallySampledMetricsTask`` for a
349 single detector.
350 bad_mask_threshold : `float`, optional
351 Samples whose summed mask fractions across ``bad_mask_columns``
352 meet or exceed this value are dropped before computing statistics
353 and rendering the figure.
354 bad_mask_columns : iterable of `str`, optional
355 Names of mask-fraction columns to sum when deciding whether a
356 sample sits on a usable patch of the detector. Columns missing
357 from the input table are silently ignored.
358 chi2_scale : `float`, optional
359 Multiplicative half-width of the ``diffim_chi2PerPix`` colormap
360 around its nominal value of 1.0. The colorbar covers
361 ``[1/(1+chi2_scale), 1+chi2_scale]`` so a deviation by a factor
362 of ``(1+chi2_scale)`` above or below 1 sits at the colormap extremes.
363 skymap : `lsst.skymap.BaseSkyMap`, optional
364 If supplied, overlay the boundaries of every patch that touches
365 the detector footprint on each panel, with a ``tract,patch``
366 label anchored just inside the lower-left corner of each patch.
367 The local sky↔pixel mapping is inferred from the sample
368 positions' ``(x, y)`` and ``(coord_ra, coord_dec)`` columns via
369 a least-squares affine fit, so the alignment is approximate
370 When omitted, no outlines are drawn.
371 label_fontsize : `int` or `float`, optional
372 Font size for the per-patch ``tract,patch`` labels. Only used
373 when ``skymap`` is supplied. Default 7.
374 panel_padding_pix : `float`, optional
375 Detector-pixel buffer added on each side of every panel beyond
376 the autoscaled data limits.
378 Returns
379 -------
380 fig : `matplotlib.figure.Figure`
381 The 1x3 diagnostic figure (`diffim_chi2PerPix`, PSF match
382 residual, dipole density). The kernel centroid offset quiver is
383 overlaid on every panel, with a quiverkey reference arrow.
385 Notes
386 -----
387 The headline summary is printed via ``print``, so the function is
388 intended for direct use in a notebook cell. Pass the returned figure
389 to ``fig.savefig(...)`` if you want to persist the diagnostic.
390 """
391 df = _coerce_to_frame(metrics)
392 clean = _filter_clean(df, bad_mask_threshold, bad_mask_columns)
393 _print_summary(clean, n_total=len(df), threshold=bad_mask_threshold)
394 return _make_figure(clean, chi2_scale=chi2_scale, skymap=skymap,
395 label_fontsize=label_fontsize,
396 panel_padding_pix=panel_padding_pix)