Coverage for python/lsst/analysis/ap/nb_utils.py: 12%
885 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-23 09:35 +0000
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/>.
22from __future__ import annotations
24__all__ = ["make_simbad_link", "compare_sources", "compare_objects",
25 "find_objects_sharing_sources",
26 "classify_association_clusters",
27 "plot_cutouts_with_object_markers",
28 "plot_objects_sharing_sources",
29 "display_images", "display_images_ab", "display_footprints",
30 "display_coadd_coverage",
31 "get_xy_from_source_table", "extract_timestamped_messages"]
33import astropy.coordinates as coord
34from astroquery.simbad import Simbad
35import astropy.units as u
36import astropy.table
37from dataclasses import dataclass, field
38from datetime import datetime, timezone
39import functools
40import json
41import numpy as np
42import os
43import random
44import pandas as pd
45from typing import Any
47import lsst.afw.display
48import lsst.afw.table
49import lsst.geom
50from lsst.daf.butler import DatasetNotFoundError
51from lsst.images import Box
52from lsst.analysis.ap import plotImageSubtractionCutouts
53from lsst.analysis.ap.compare import match_catalogs
54from lsst.analysis.ap.imageUtils import (
55 _box_edges, _pixel_to_sky, _sky_to_pixel, cutout, get_visit_image,
56)
57from lsst.analysis.ap.skymapOverlay import draw_skymap_outlines_afw
58# Polygon geometry shared with the skymap overlays; private to the
59# package, not to the module.
60from lsst.analysis.ap.skymapOverlay import _clip_polygon_to_rect, _polygon_area
61from IPython.display import display, Image, Markdown
64# Maps the image_type kwarg used by `display_images` to butler dataset
65# names.
66_IMAGE_DATASETS = {
67 "science": "preliminary_visit_image",
68 "template": "template_detector",
69 "difference": "difference_image",
70}
73def _apply_fakes_prefix(image_datasets, use_fakes):
74 """Prepend the fake-source pipeline's prefixes to the science and
75 template image dataset names when ``use_fakes`` is True: ``fakes_``
76 for the science image, ``injectedTemplate_`` for the template. The
77 difference image and all catalogs keep their non-prefixed names
78 (the fake-source pipeline injects into science + template but
79 re-uses the same downstream artifacts).
80 """
81 if not use_fakes:
82 return image_datasets
83 return {
84 "science": f"fakes_{image_datasets['science']}",
85 "template": f"injectedTemplate_{image_datasets['template']}",
86 "difference": image_datasets["difference"],
87 }
90def _cutout_exists(cpath, dia_source_id):
91 """Return True if a cutout PNG for this diaSourceId already exists.
93 Parameters
94 ----------
95 cpath : `~plotImageSubtractionCutouts.CutoutPath`
96 Path manager for the cutout directory.
97 dia_source_id : `int`
98 DiaSourceId whose ``{id}.png`` is checked.
99 """
100 return cpath.exists(dia_source_id, f"{dia_source_id}.png")
103def make_simbad_link(ra, dec, radius_arcsec=3.0):
104 """Search Simbad for associated sources within a 3 arcsecond region.
106 Parameters
107 ----------
108 ra : 'float'
109 Ra from source.
111 dec : 'float'
112 Dec from source.
114 radius_arcsec : 'float'
115 Search radius submitted to Simbad in arcseconds.
116 Default radius is 3 arcseconds.
118 Returns
119 -------
120 results_table : `astropy.table.table.Table`
121 A table of Simbad search results.
122 """
123 search_results = f"http://simbad.cds.unistra.fr/simbad/sim-coo?Coord={ra}+{dec}" \
124 f"&CooFrame=FK5&CooEpoch=2000&CooEqui=2000&CooDefinedFrames=none&Radius=" \
125 f"{radius_arcsec}&Radius.unit=arcsec&submit=submit+query&CoordList="
126 display(Markdown(f"[Link to Simbad search]({search_results})"))
128 source_coords = coord.SkyCoord(ra, dec, frame="icrs", unit=(u.deg, u.deg))
129 customSimbad = Simbad()
130 customSimbad.TIMEOUT = 600
131 customSimbad.add_votable_fields("otype(V)")
132 results_table = customSimbad.query_region(
133 source_coords, radius=radius_arcsec*u.arcsecond
134 )
136 if results_table is not None:
138 return results_table
140 else:
141 print(f"No matched sources within {radius_arcsec} arcseconds.")
143 return None
146def compare_sources(butler1, butler2, query1, query2,
147 bad_flag_list=None, match_radius=0.5,
148 make_cutouts=False, display_cutouts=False,
149 cutout_path1=None, cutout_path2=None,
150 cutout_config1=None, cutout_config2=None,
151 njobs=0):
152 """Compare two APDB datasets by extracting unassociated sources,
153 spatially crossmatching, and plotting cutouts of the differences.
155 Parameters
156 ----------
157 butler1 : `lsst.daf.butler`
158 Initialized Butler repo containing the first dataset.
159 Could be the same as butler2 but should be initialized with the
160 appropriate collection name for cutout generation, if doing that.
161 butler2 : `lsst.daf.butler`
162 Initialized Butler repo containing the second dataset.
163 Could be the same as butler1 but should be initialized with the
164 appropriate collection name for cutout generation, if doing that.
165 query1 : `lsst.analysis.ap.DbQuery`
166 DbQuery to first APDB (postgresql or slite file;
167 NOT created in this function).
168 query2 : `lsst.analysis.ap.DbQuery`
169 DbQuery to second APDB (postgresql or slite file;
170 NOT created in this function).
171 bad_flag_list : `list`, optional
172 List of bad flags to exclude (applied to both query1 and query2).
173 Omit list to skip filtering.
174 match_radius : `double`
175 Maximum allowable distance in arcsec between an object in
176 data1 and data2.
177 make_cutouts : `bool`, optional
178 Generate cutouts for sources unique to each dataset; default is False.
179 display_cutouts: `bool`, optional
180 Display cutouts for sources present in only one of the DBs to the
181 screen; default is False.
182 cutout_path1, cutout_path2 : `str`, optional
183 Base path to store cutouts for sources unique to the datasets.
184 Must be supplied if make_cutouts is True.
185 cutout_config1, cutout_config2 : `dict` [`str`], optional
186 Config overrides to apply to cutout plotter for the datasets.
187 See `~plotImageSubtractionCutouts.PlotImageSubtractionCutoutsConfig`
188 for available options.
189 njobs : `int`, optional
190 Number of parallel processes for plotImageSubtractionCutouts.
192 Returns
193 -------
194 unique1 : `pandas.DataFrame`
195 Data frame of sources only found in the first dataset.
196 unique2 : `pandas.DataFrame`
197 Data frame of sources only found in the second dataset.
198 matched : `pandas.DataFrame`
199 Data frame of matched sources; the rows are sources from the first
200 dataset, with two columns added: ``src2_diaSourceId`` pointing to
201 the matched diaSourceId in the second dataset, and
202 ``xmatch_dist_arcsec`` giving the on-sky separation in arcseconds.
203 """
205 if make_cutouts and (cutout_path1 is None or cutout_path2 is None):
206 errstr = ('You must supply a value for `cutout_path1` and `cutout_path2` if `make_cutouts` is True.')
207 raise ValueError(errstr)
209 if bad_flag_list is not None:
210 # Snapshot and restore so we don't leave the caller's queries with a
211 # different exclusion list than they started with.
212 saved_flags1 = list(query1.diaSource_flags_exclude)
213 saved_flags2 = list(query2.diaSource_flags_exclude)
214 query1.set_excluded_diaSource_flags(bad_flag_list)
215 query2.set_excluded_diaSource_flags(bad_flag_list)
216 try:
217 goodSrc1 = query1.load_sources(exclude_flagged=True)
218 goodSrc2 = query2.load_sources(exclude_flagged=True)
219 finally:
220 query1.set_excluded_diaSource_flags(saved_flags1)
221 query2.set_excluded_diaSource_flags(saved_flags2)
222 else:
223 goodSrc1 = query1.load_sources(exclude_flagged=True)
224 goodSrc2 = query2.load_sources(exclude_flagged=True)
226 if 'reliability' not in goodSrc1.columns:
227 goodSrc1['reliability'] = None
228 if 'reliability' not in goodSrc2.columns:
229 goodSrc2['reliability'] = None
231 # Cross-match within each (visit, detector) group.
232 matched, unique1, unique2 = match_catalogs(
233 goodSrc1, goodSrc2,
234 radius=match_radius * u.arcsec,
235 on=("visit", "detector"),
236 )
237 # Preserve the legacy column name on the returned `matched` DataFrame.
238 matched = matched.rename(columns={"diaSourceId_2": "src2_diaSourceId"})
240 print("{} matched sources; {} unique to set 1; {} unique to set 2.".format(
241 len(matched), len(unique1), len(unique2)))
243 # Decide if we are doing anything with cutouts or not. If not, just skip.
244 if make_cutouts:
245 # Make paths if they don't exist.
246 if not os.path.exists(cutout_path1):
247 os.makedirs(cutout_path1)
248 if not os.path.exists(cutout_path2):
249 os.makedirs(cutout_path2)
251 # Make cutouts if they don't already exist
252 config1 = plotImageSubtractionCutouts.PlotImageSubtractionCutoutsConfig()
253 config2 = plotImageSubtractionCutouts.PlotImageSubtractionCutoutsConfig()
254 # default to flat directories for ease of use
255 config1.chunk_size = None
256 config2.chunk_size = None
257 # apply user-specified overrides
258 if cutout_config1 is not None:
259 config1.update(**cutout_config1)
260 if cutout_config2 is not None:
261 config2.update(**cutout_config2)
263 cpath1 = plotImageSubtractionCutouts.CutoutPath(cutout_path1,
264 chunk_size=config1.chunk_size)
265 cpath2 = plotImageSubtractionCutouts.CutoutPath(cutout_path2,
266 chunk_size=config2.chunk_size)
268 plotter1 = plotImageSubtractionCutouts.PlotImageSubtractionCutoutsTask(
269 output_path=cutout_path1, config=config1)
270 plotter2 = plotImageSubtractionCutouts.PlotImageSubtractionCutoutsTask(
271 output_path=cutout_path2, config=config2)
273 # First figure out which cutouts already exist at the output path.
274 # Series.apply passes one positional argument (the diaSourceId), but
275 # _cutout_exists also needs the per-dataset cpath; partial binds it.
276 unique1['pathexists'] = unique1['diaSourceId'].apply(
277 functools.partial(_cutout_exists, cpath1))
278 pathchk1 = unique1.loc[~unique1['pathexists']]
280 unique2['pathexists'] = unique2['diaSourceId'].apply(
281 functools.partial(_cutout_exists, cpath2))
282 pathchk2 = unique2.loc[~unique2['pathexists']]
284 # Only write those that don't exist yet
285 plotter1.write_images(pathchk1, butler=butler1, njobs=njobs)
286 plotter2.write_images(pathchk2, butler=butler2, njobs=njobs)
288 if display_cutouts:
289 for isrc in unique1.itertuples():
290 fpath = cpath1(int(isrc.diaSourceId), f"{int(isrc.diaSourceId)}.png")
291 print('Unique to dataset 1: {}'.format(int(isrc.diaSourceId)))
292 display(Image(filename=fpath))
294 for isrc in unique2.itertuples():
295 fpath = cpath2(int(isrc.diaSourceId), f"{int(isrc.diaSourceId)}.png")
297 print('Unique to dataset 2: {}'.format(int(isrc.diaSourceId)))
298 display(Image(filename=fpath))
300 # drop pathexists columns to return to original dataframe shape
301 _ = unique1.pop('pathexists')
302 _ = unique2.pop('pathexists')
304 return unique1, unique2, matched
307def compare_objects(query1, query2, match_radius=0.5):
308 """Compare two APDB datasets by spatially crossmatching diaObjects.
310 Parameters
311 ----------
312 query1 : `lsst.analysis.ap.DbQuery`
313 DbQuery to first APDB (postgresql or slite file;
314 NOT created in this function).
315 query2 : `lsst.analysis.ap.DbQuery`
316 DbQuery to second APDB (postgresql or slite file;
317 NOT created in this function).
318 match_radius : `double`
319 Maximum allowable distance in arcsec between an object in
320 data1 and data2.
322 Returns
323 -------
324 unique1 : `pandas.DataFrame`
325 Data frame of diaObjects only found in the first dataset.
326 unique2 : `pandas.DataFrame`
327 Data frame of diaObjects only found in the second dataset.
328 matched : `pandas.DataFrame`
329 Data frame of matched diaObjects; the rows are objects from the
330 first dataset, with two columns added: ``obj2_diaObjectId``
331 pointing to the matched diaObjectId in the second dataset, and
332 ``xmatch_dist_arcsec`` giving the on-sky separation in arcseconds.
333 """
334 obj1 = query1.load_objects()
335 obj2 = query2.load_objects()
337 # diaObjects aren't tied to a single (visit, detector); match across
338 # the full catalog with `on=()`.
339 matched, unique1, unique2 = match_catalogs(
340 obj1, obj2,
341 radius=match_radius * u.arcsec,
342 on=(),
343 id_col="diaObjectId",
344 )
345 matched = matched.rename(columns={"diaObjectId_2": "obj2_diaObjectId"})
347 print("{} matched objects; {} unique to set 1; {} unique to set 2.".format(
348 len(matched), len(unique1), len(unique2)))
350 return unique1, unique2, matched
353def _match_source_ids(sources1, sources2, match_radius):
354 """Return the (diaSourceId, diaSourceId_2) correspondence between two
355 runs' diaSource catalogs, plus run 1's sky position for each pair.
357 diaSourceIds are not stable across runs -- they end in a per-catalog
358 counter assigned in detection order, which any change to detection
359 or measurement renumbers -- so the same detection must be identified
360 by position within a single (visit, detector).
362 Parameters
363 ----------
364 sources1, sources2 : `pandas.DataFrame`
365 diaSource catalogs, each with `diaSourceId`, `diaObjectId`,
366 `ra`, `dec`, `visit` and `detector`.
367 match_radius : `float`
368 Maximum separation in arcsec for a pair to count as the same
369 detection.
371 Returns
372 -------
373 paired : `pandas.DataFrame`
374 Columns `diaSourceId`, `diaObjectId`, `ra`, `dec` (all from run
375 1), `diaSourceId_2` and `diaObjectId_2` (from run 2). One row
376 per matched pair; unmatched sources are absent.
377 """
378 cols = ["diaSourceId", "diaObjectId", "ra", "dec", "visit", "detector"]
379 matched, _, _ = match_catalogs(sources1[cols], sources2[cols],
380 radius=match_radius * u.arcsec,
381 on=("visit", "detector"))
382 # match_catalogs gives every run-1 source its nearest run-2 neighbor, so
383 # two run-1 sources can claim the same run-2 source; keep only the closest
384 # of those so the pairing stays one-to-one.
385 matched = matched.sort_values("xmatch_dist_arcsec").drop_duplicates(
386 subset="diaSourceId_2")
387 return matched[["diaSourceId", "diaObjectId", "ra", "dec",
388 "diaSourceId_2"]].merge(
389 sources2[["diaSourceId", "diaObjectId"]].rename(
390 columns={"diaSourceId": "diaSourceId_2",
391 "diaObjectId": "diaObjectId_2"}),
392 on="diaSourceId_2", how="inner")
395def _to_run1_ids(sources2, obj2_ids, id2_to_id1):
396 """Return the run-1 diaSourceIds of every run-2 diaSource owned by
397 one of ``obj2_ids``, dropping any with no run-1 counterpart.
398 """
399 ids2 = sources2.loc[sources2["diaObjectId"].isin(obj2_ids), "diaSourceId"]
400 return {id2_to_id1[i] for i in ids2 if i in id2_to_id1}
403def find_objects_sharing_sources(diaObjectId, sources1, sources2,
404 objects1, objects2,
405 max_distance_arcsec=2, match_radius=0.5):
406 """For a diaObjectId in run 1, return the full association cluster
407 of diaSources and diaObjects from both runs.
409 Treats the (diaSource, run-1-diaObject, run-2-diaObject) links as a
410 graph -- each diaSource is connected to its owning diaObject in
411 each run -- and grows the connected component reachable from the
412 input diaObjectId until no new sources or objects are discovered.
413 Catches arbitrarily deep merge/split chains across the two runs
414 (e.g. run 2 merges A+B into Z, then a third source in B is split
415 into a fourth object in run 2, etc.).
417 diaSourceIds are *not* stable across runs (they end in a per-catalog
418 counter assigned in detection order), so the two runs' diaSources
419 are paired by position via `_match_source_ids`. Only detections
420 present in both runs carry the graph; a diaSource with no
421 counterpart within ``match_radius`` cannot link objects across runs.
423 Parameters
424 ----------
425 diaObjectId : `int`
426 A diaObjectId, typically from `unique1` returned by
427 `compare_objects`.
428 sources1, sources2 : `pandas.DataFrame`
429 Full diaSources catalogs from runs 1 and 2 (e.g. from
430 ``query.load_sources()``). Each must contain `diaSourceId`,
431 `diaObjectId`, `ra`, and `dec` columns.
432 objects1, objects2 : `pandas.DataFrame`
433 Full diaObjects catalogs from runs 1 and 2 (e.g. from
434 ``query.load_objects()``). Each must contain `diaObjectId`,
435 `ra`, and `dec` columns.
436 max_distance_arcsec : `float`, optional
437 If given, only include diaSources within this distance of the input
438 diaObject's (ra, dec) in the search. All diaSources of the final
439 diaObjects will still be returned, even if outside this distance.
440 match_radius : `float`, optional
441 Maximum separation in arcsec for a run-1 and a run-2 diaSource to
442 be treated as the same detection.
444 Returns
445 -------
446 sources : `pandas.DataFrame`
447 Rows of `sources1` for every diaSource belonging to any of the
448 found run-2 diaObjects (in run 2's view).
449 related_objects1 : `pandas.DataFrame`
450 Rows of `objects1` for every run-1 diaObject containing any of
451 those diaSources in run 1.
452 related_objects2 : `pandas.DataFrame`
453 Rows of `objects2` for every run-2 diaObject containing any of
454 those diaSources in run 2.
455 """
456 # Pair the two runs' diaSources up front; the search below runs in
457 # run-1 id space and translates run-2 sources through this map.
458 paired = _match_source_ids(sources1, sources2, match_radius)
459 id2_to_id1 = dict(zip(paired["diaSourceId_2"], paired["diaSourceId"]))
460 id1_to_id2 = dict(zip(paired["diaSourceId"], paired["diaSourceId_2"]))
462 if max_distance_arcsec is not None:
463 ref_match = objects1[objects1["diaObjectId"] == diaObjectId]
464 if len(ref_match) == 0:
465 raise ValueError(
466 f"diaObjectId={diaObjectId} not found in objects1")
467 ref_row = ref_match.iloc[0]
468 ref = coord.SkyCoord(ra=ref_row["ra"] * u.deg,
469 dec=ref_row["dec"] * u.deg)
470 # The search space is run-1 diaSourceIds, so filter against
471 # sources1; positions agree between paired sources to within
472 # match_radius.
473 sep = ref.separation(
474 coord.SkyCoord(ra=sources1["ra"].values * u.deg,
475 dec=sources1["dec"].values * u.deg)
476 ).to_value(u.arcsec)
477 allowed_src_ids = set(
478 sources1.loc[sep <= max_distance_arcsec, "diaSourceId"])
479 else:
480 allowed_src_ids = None
482 # Breadth-first search for the connected component:
483 # alternately expand sources from the currently-known objects,
484 # then expand objects from the sources.
485 # Terminates because every iteration adds at least one source
486 # before the fixed-point check fires, and the source pool is finite.
487 src_ids = set()
488 obj1_ids = {diaObjectId}
489 obj2_ids = set()
491 while True:
492 new_src_ids = set(
493 sources1.loc[sources1["diaObjectId"].isin(obj1_ids),
494 "diaSourceId"])
495 new_src_ids.update(_to_run1_ids(sources2, obj2_ids, id2_to_id1))
496 if allowed_src_ids is not None:
497 new_src_ids &= allowed_src_ids
498 if new_src_ids <= src_ids:
499 break
500 src_ids |= new_src_ids
501 obj1_ids |= set(
502 sources1.loc[sources1["diaSourceId"].isin(src_ids),
503 "diaObjectId"])
504 run2_ids = {id1_to_id2[i] for i in src_ids if i in id1_to_id2}
505 obj2_ids |= set(
506 sources2.loc[sources2["diaSourceId"].isin(run2_ids),
507 "diaObjectId"])
509 # Expand the final diaSource list to every source owned by any
510 # surviving diaObject.
511 final_src_ids = set(
512 sources1.loc[sources1["diaObjectId"].isin(obj1_ids), "diaSourceId"])
513 final_src_ids |= _to_run1_ids(sources2, obj2_ids, id2_to_id1)
515 sources = sources1[sources1["diaSourceId"].isin(final_src_ids)]
516 related_objects1 = objects1[objects1["diaObjectId"].isin(obj1_ids)]
517 related_objects2 = objects2[objects2["diaObjectId"].isin(obj2_ids)]
519 return sources, related_objects1, related_objects2
522class _UnionFind:
523 """Disjoint-set with path compression and union-by-rank.
525 Used by `classify_association_clusters` to quickly find connected
526 components of the (run-1 diaObject, run-2 diaObject) graph.
527 """
529 def __init__(self):
530 self._parent = {}
531 self._rank = {}
533 def add(self, x):
534 if x not in self._parent:
535 self._parent[x] = x
536 self._rank[x] = 0
538 def find(self, x):
539 # Two-pass iterative find with path compression.
540 root = x
541 while self._parent[root] != root:
542 root = self._parent[root]
543 while self._parent[x] != root:
544 self._parent[x], x = root, self._parent[x]
545 return root
547 def union(self, x, y):
548 rx, ry = self.find(x), self.find(y)
549 if rx == ry:
550 return
551 if self._rank[rx] < self._rank[ry]:
552 rx, ry = ry, rx
553 self._parent[ry] = rx
554 if self._rank[rx] == self._rank[ry]:
555 self._rank[rx] += 1
558def classify_association_clusters(sources1, sources2, match_radius=0.5):
559 """Enumerate and classify every association-disagreement cluster
560 between two APDBs that share input diaSources.
562 Builds the bipartite graph whose edges are
563 ``(diaSource -> its run-1 diaObject, diaSource -> its run-2
564 diaObject)`` over all diaSources the two runs have in common, runs
565 union-find over the diaObjectIds to extract every connected
566 component, and labels each cluster:
568 * ``matched`` -- one run-1 obj <-> one run-2 obj.
569 * ``split`` -- one run-1 obj split into multiple run-2 objs.
570 * ``merged`` -- multiple run-1 objs merged into one run-2 obj.
571 * ``tangled`` -- M run-1 objs <-> N run-2 objs, both > 1.
573 diaSourceIds are *not* stable across runs: they carry a per-catalog
574 counter assigned in detection order, so any change to detection or
575 measurement renumbers them. The common diaSources are therefore
576 identified by position (nearest neighbor within ``match_radius``,
577 inside a single (visit, detector)) rather than by id. Sources with
578 no counterpart in the other run are skipped.
580 Parameters
581 ----------
582 sources1, sources2 : `pandas.DataFrame`
583 Full diaSources catalogs from runs 1 and 2 (e.g. from
584 ``query.load_sources()``). Each must contain `diaSourceId`,
585 `diaObjectId`, `ra`, `dec`, `visit`, and `detector` columns.
586 match_radius : `float`, optional
587 Maximum separation in arcsec for two diaSources to be considered
588 the same detection in both runs.
590 Returns
591 -------
592 clusters : `pandas.DataFrame`
593 One row per cluster, with columns:
594 - ``kind``: matched / split / merged / tangled.
595 - ``n_obj1``, ``n_obj2``: distinct diaObject counts per run.
596 - ``n_sources``: matched diaSource pairs in the cluster.
597 - ``obj1_ids``, ``obj2_ids``: tuples of diaObjectIds.
598 - ``ra``, ``dec``: mean sky position of the cluster's
599 diaSources (degrees).
600 """
601 # Pre-define the types so that value_counts() and groupby()
602 # include unused kinds with a count of 0.
603 kind_dtype = pd.CategoricalDtype(
604 categories=["matched", "split", "merged", "tangled"], ordered=True)
606 paired = _match_source_ids(sources1, sources2, match_radius)
608 if len(paired) == 0:
609 empty = pd.DataFrame(columns=[
610 "kind", "n_obj1", "n_obj2", "n_sources",
611 "obj1_ids", "obj2_ids", "ra", "dec"])
612 empty["kind"] = empty["kind"].astype(kind_dtype)
613 return empty
615 # Define a namespace for the two runs so identical numeric ids in run 1 and
616 # run 2 don't collide as keys.
617 keys1 = [("r1", int(i)) for i in paired["diaObjectId"].to_numpy()]
618 keys2 = [("r2", int(i)) for i in paired["diaObjectId_2"].to_numpy()]
620 uf = _UnionFind()
621 for k in set(keys1):
622 uf.add(k)
623 for k in set(keys2):
624 uf.add(k)
625 for k1, k2 in zip(keys1, keys2):
626 uf.union(k1, k2)
628 paired = paired.assign(_cluster=[uf.find(k) for k in keys1])
630 rows = []
631 for _, grp in paired.groupby("_cluster", sort=False):
632 ids_a = tuple(sorted(int(i) for i in grp["diaObjectId"].unique()))
633 ids_b = tuple(sorted(int(i) for i in grp["diaObjectId_2"].unique()))
634 n1, n2 = len(ids_a), len(ids_b)
635 if n1 == 1 and n2 == 1:
636 kind = "matched"
637 elif n1 == 1:
638 kind = "split"
639 elif n2 == 1:
640 kind = "merged"
641 else:
642 kind = "tangled"
643 rows.append({
644 "kind": kind,
645 "n_obj1": n1, "n_obj2": n2,
646 "n_sources": grp["diaSourceId"].nunique(),
647 "obj1_ids": ids_a, "obj2_ids": ids_b,
648 "ra": float(grp["ra"].mean()),
649 "dec": float(grp["dec"].mean()),
650 })
652 result = pd.DataFrame(rows)
653 result["kind"] = result["kind"].astype(kind_dtype)
654 return result
657# Colors used by the cutout plotters to give each distinct
658# diaObjectId its own marker color.
659_OBJECT_PALETTE = ("lime", "red", "cyan", "magenta", "yellow", "orange",
660 "deepskyblue", "pink", "white", "violet", "gold",
661 "lightgreen")
664def _prepare_object_overlays(objects, palette):
665 """Deduplicate `objects` by diaObjectId and assign one palette color
666 per distinct id, returning the parallel arrays the cutout renderer
667 needs: ``(ids, ras, decs, colors)``. Done once per call so the same
668 color identifies the same diaObject across every cutout.
669 """
670 obj_unique = objects.drop_duplicates(subset="diaObjectId")
671 # Prefer the run-2 id when present (matched rows carry both); pandas
672 # concat promotes the column to float64 if any rows lack it, so cast
673 # back to int64 after filling.
674 if "obj2_diaObjectId" in obj_unique.columns:
675 obj_ids = obj_unique["obj2_diaObjectId"].combine_first(
676 obj_unique["diaObjectId"]).astype(np.int64).to_numpy()
677 else:
678 obj_ids = obj_unique["diaObjectId"].astype(np.int64).to_numpy()
679 obj_ras = np.asarray(obj_unique["ra"])
680 obj_decs = np.asarray(obj_unique["dec"])
681 obj_colors = [palette[i % len(palette)] for i in range(len(obj_ids))]
682 return obj_ids, obj_ras, obj_decs, obj_colors
685def _load_cutout(butler, row, *, size, image_type, image_datasets):
686 """Fetch the requested image dataset for this row's (visit,
687 detector) and return a small dict with everything the renderer
688 needs: pixel data, dimensions, cutout origin, WCS, an
689 ImageNormalize tuned to the central source, and the `image_type`
690 label used in the cutout title.
692 Loading is separated from rendering so callers that need to draw
693 the same cutout into multiple Axes can pay the read and cutout cost
694 once.
695 """
696 import astropy.visualization as aviz
698 data_id = {"visit": int(row.visit), "detector": int(row.detector)}
699 image = _get_image(butler, image_datasets, image_type, data_id)
700 stamp = cutout(image, coord.SkyCoord(row.ra, row.dec, unit="deg"), size)
701 data = stamp.array
702 ny, nx = data.shape
704 if image_type == "difference": 704 ↛ 711line 704 didn't jump to line 711 because the condition on line 704 was always true
705 # Normalize on a small central window so the source dominates
706 # the dynamic range.
707 cy, cx = ny // 2, nx // 2
708 half = min(7, cy, cx)
709 norm_data = data[cy - half:cy + half + 1, cx - half:cx + half + 1]
710 else:
711 norm_data = data
712 norm = aviz.ImageNormalize(
713 norm_data, interval=aviz.MinMaxInterval(),
714 stretch=aviz.AsinhStretch(a=0.1))
716 return {
717 "data": data, "ny": ny, "nx": nx,
718 "sky_projection": image.sky_projection,
719 "x0": stamp.bbox.x.start, "y0": stamp.bbox.y.start,
720 "norm": norm,
721 "image_type": image_type,
722 }
725def _render_cutout_axes(ax, row, cutout_data, sources,
726 obj_ids, obj_ras, obj_decs, obj_colors, *,
727 marker_size, marker_symbol,
728 source_marker_size, current_source_marker_size,
729 current_source_color,
730 source_match_ids=None,
731 title=None, subtitle=""):
732 """Render one diaSource cutout onto an existing matplotlib Axes
733 using preloaded data from `_load_cutout`.
735 Internal helper shared by `plot_cutouts_with_object_markers` and
736 `plot_objects_sharing_sources`. The caller owns figure creation,
737 layout, saving, and displaying.
739 By default the axes title is built from `row` as
740 ``"diaSourceId=... (image_type, visit=..., det=...)"``. Pass
741 `title=` explicitly (including ``""``) to override or suppress
742 that line -- useful when a parent figure or subfigure already
743 carries the shared header. If `subtitle` is non-empty it is drawn
744 on a second title line.
745 """
746 from matplotlib import cm
748 data = cutout_data["data"]
749 ny = cutout_data["ny"]
750 nx = cutout_data["nx"]
751 x0 = cutout_data["x0"]
752 y0 = cutout_data["y0"]
753 sky_projection = cutout_data["sky_projection"]
754 norm = cutout_data["norm"]
755 image_type = cutout_data["image_type"]
757 ax.imshow(data, cmap=cm.bone, interpolation="none", norm=norm,
758 origin="lower", aspect="equal",
759 extent=(0, nx, 0, ny))
761 this_id = int(row.diaSourceId)
763 # Project every supplied diaSource into the cutout frame once,
764 # then split into "this cutout's diaSource" vs every other
765 # diaSource whose sky position falls inside this cutout
766 # (regardless of which image it was detected on).
767 if len(sources) > 0:
768 src_xs, src_ys = _sky_to_pixel(sky_projection, sources["ra"], sources["dec"])
769 src_xs = src_xs - x0
770 src_ys = src_ys - y0
771 src_in_frame = (
772 (src_xs >= 0) & (src_xs < nx)
773 & (src_ys >= 0) & (src_ys < ny))
774 id_arr = sources["diaSourceId"].to_numpy()
775 other_src_mask = src_in_frame & (id_arr != this_id)
776 current_src_mask = src_in_frame & (id_arr == this_id)
777 else:
778 src_xs = src_ys = None
779 other_src_mask = current_src_mask = None
781 # Per-diaSource color resolution: map each diaSource to its
782 # owning diaObject (in this panel's view), then to that
783 # diaObject's palette color. Falls back to `current_source_color`
784 # for sources whose owner is not present in `obj_ids`.
785 if source_match_ids is None:
786 match_ids_list = sources["diaObjectId"].tolist()
787 else:
788 match_ids_list = list(source_match_ids)
789 src_to_match = dict(zip(sources["diaSourceId"].tolist(),
790 match_ids_list))
791 id_to_color = dict(zip(obj_ids.tolist(), obj_colors))
793 def _color_for(diaSourceId):
794 return id_to_color.get(
795 src_to_match.get(int(diaSourceId)), current_source_color)
797 if other_src_mask is not None and other_src_mask.any():
798 other_indices = np.flatnonzero(other_src_mask)
799 other_colors = [_color_for(int(id_arr[i])) for i in other_indices]
800 ax.scatter(src_xs[other_src_mask], src_ys[other_src_mask],
801 s=source_marker_size, marker="+",
802 c=other_colors, linewidths=1.0,
803 label="other diaSource")
805 if len(obj_ids) > 0:
806 xs, ys = _sky_to_pixel(sky_projection, obj_ras, obj_decs)
807 xs = xs - x0
808 ys = ys - y0
809 in_bounds = (xs >= 0) & (xs < nx) & (ys >= 0) & (ys < ny)
810 else:
811 in_bounds = np.zeros(0, dtype=bool)
813 for i in np.flatnonzero(in_bounds):
814 ax.scatter(xs[i], ys[i],
815 s=marker_size, marker=marker_symbol,
816 facecolors="none", edgecolors=obj_colors[i],
817 linewidths=1.5,
818 label=f"diaObjectId={int(obj_ids[i])}")
820 # Current diaSource last so it stays on top of any overlapping
821 # diaObject marker at the cutout center.
822 if current_src_mask is not None and current_src_mask.any():
823 ax.scatter(src_xs[current_src_mask], src_ys[current_src_mask],
824 s=current_source_marker_size, marker="x",
825 c=_color_for(this_id), linewidths=2.0,
826 label=f"current diaSourceId={this_id}")
828 if title is None:
829 title = (f"diaSourceId={this_id} "
830 f"({image_type}, visit={int(row.visit)}, "
831 f"det={int(row.detector)})")
832 if subtitle:
833 title = f"{title}\n{subtitle}" if title else subtitle
834 if title:
835 ax.set_title(title, fontsize="small")
836 ax.set_xticks([])
837 ax.set_yticks([])
838 if ax.get_legend_handles_labels()[0]:
839 ax.legend(loc="upper right", fontsize="x-small", framealpha=0.7)
842def plot_cutouts_with_object_markers(sources, butler, objects, *,
843 output_path=None,
844 display_cutouts=False,
845 size=51,
846 image_type="difference",
847 image_datasets=_IMAGE_DATASETS,
848 marker_size=80,
849 marker_symbol="o",
850 palette=_OBJECT_PALETTE,
851 source_marker_size=80,
852 current_source_marker_size=180,
853 current_source_color="yellow"):
854 """Plot per-diaSource cutouts with overlaid markers at given diaObject
855 sky positions.
857 For each diaSource in `sources`, fetch a square cutout from `butler`
858 centered on the source's (ra, dec). On each cutout draw:
860 * A small ``+`` marker at every other diaSource in `sources`
861 whose sky position lands inside the cutout, regardless of which
862 (visit, detector) it was detected on.
863 * A distinct ``x`` marker for the diaSource the cutout is
864 centered on (the "current" diaSource).
865 * One color-coded marker per distinct diaObjectId in `objects`,
866 cycling through `palette`; the same color identifies the same
867 diaObject across every cutout in the run.
869 Markers that fall outside the cutout bounds are skipped.
871 Typical use: visualize how a group of diaSources (all originally
872 associated with one diaObjectId in run 1) got redistributed across
873 diaObjects in run 2. `sources` and `objects` are usually built from
874 the output of `find_objects_sharing_sources`::
876 sources, ro1, ro2 = find_objects_sharing_sources(
877 diaObjectId, sources1, sources2, objects1, objects2)
878 objects = pd.concat([ro1, ro2])
879 plot_cutouts_with_object_markers(
880 sources, butler1, objects, display_cutouts=True,
881 )
883 Parameters
884 ----------
885 sources : `pandas.DataFrame`
886 DiaSources to cut out. Must contain `diaSourceId`, `ra`, `dec`,
887 `visit`, and `detector` columns.
888 butler : `lsst.daf.butler.Butler`
889 Butler containing the image datasets for these (visit, detector)
890 pairs.
891 objects : `pandas.DataFrame`
892 DiaObjects to mark. Must contain `diaObjectId`, `ra`, and `dec`
893 columns. Duplicate diaObjectIds are dropped (first row wins).
894 If an `obj2_diaObjectId` column is present (e.g. for rows from
895 a `matched` DataFrame returned by `compare_objects`), the run-2
896 id is shown in the legend in preference to the run-1
897 `diaObjectId`.
898 output_path : `str`, optional
899 Directory to write ``{diaSourceId}.png`` files to. Created if
900 missing. Pass None to skip writing.
901 display_cutouts : `bool`, optional
902 If True, display each cutout inline (notebook).
903 size : `int`, optional
904 Cutout side length in pixels.
905 image_type : {"science", "template", "difference"}, optional
906 Which image to render.
907 image_datasets : `dict` [`str`, `str`], optional
908 Mapping from image-type key to butler dataset name.
909 marker_size : `int`, optional
910 matplotlib scatter ``s`` parameter for diaObject markers.
911 marker_symbol : `str`, optional
912 matplotlib scatter ``marker`` parameter for diaObject markers.
913 palette : sequence of `str`, optional
914 Color cycle used to assign one color per diaObjectId.
915 source_marker_size : `int`, optional
916 Scatter ``s`` parameter for the small ``+`` markers drawn at
917 the positions of the *other* diaSources in `sources`.
918 current_source_marker_size : `int`, optional
919 Scatter ``s`` parameter for the distinct marker drawn at the
920 diaSource the cutout is centered on.
921 current_source_color : `str`, optional
922 Color of the current-diaSource marker.
923 """
924 import matplotlib.pyplot as plt
926 if image_type not in image_datasets:
927 raise ValueError(
928 f"image_type must be one of {sorted(image_datasets)}, "
929 f"got {image_type!r}")
931 if output_path is not None:
932 os.makedirs(output_path, exist_ok=True)
934 overlays = _prepare_object_overlays(objects, palette)
936 for row in sources.itertuples(index=False):
937 cutout_data = _load_cutout(butler, row, size=size,
938 image_type=image_type,
939 image_datasets=image_datasets)
940 fig, ax = plt.subplots()
941 _render_cutout_axes(
942 ax, row, cutout_data, sources, *overlays,
943 marker_size=marker_size, marker_symbol=marker_symbol,
944 source_marker_size=source_marker_size,
945 current_source_marker_size=current_source_marker_size,
946 current_source_color=current_source_color)
948 if output_path is not None:
949 fpath = os.path.join(output_path, f"{int(row.diaSourceId)}.png")
950 fig.savefig(fpath, bbox_inches="tight")
951 if display_cutouts:
952 display(fig)
953 plt.close(fig)
956def plot_objects_sharing_sources(diaObjectId, sources1, sources2,
957 objects1, objects2, butler, *,
958 max_distance_arcsec=None,
959 output_path=None,
960 display_figure=True,
961 column_labels=("run 1", "run 2"),
962 figsize_per_row=4.0,
963 size=51,
964 image_type="difference",
965 image_datasets=_IMAGE_DATASETS,
966 marker_size=80,
967 marker_symbol="o",
968 palette=_OBJECT_PALETTE,
969 source_marker_size=80,
970 current_source_marker_size=180,
971 current_source_color="yellow"):
972 """Two-column cutout figure comparing the run-1 and run-2 views of
973 an association cluster.
975 Calls `find_objects_sharing_sources` internally to identify the
976 cluster of diaSources and diaObjects reachable from the input
977 `diaObjectId`, then renders one row per diaSource in the cluster.
978 The same cutout image is loaded once per row and drawn into both
979 columns: the left panel is overlaid with run-1 diaObject markers
980 (from `objects1`) and the right panel with run-2 diaObject markers
981 (from `objects2`). Each column's diaObjects get their own palette
982 mapping, so the same color in the left and right columns does
983 *not* imply the same diaObject.
985 Parameters
986 ----------
987 diaObjectId : `int`
988 Starting diaObjectId for the cluster walk.
989 sources1, sources2, objects1, objects2 : `pandas.DataFrame`
990 Forwarded to `find_objects_sharing_sources`.
991 butler : `lsst.daf.butler.Butler`
992 Butler used to fetch the cutout images. The same image backs
993 both panels of a given row.
994 max_distance_arcsec : `float`, optional
995 Forwarded to `find_objects_sharing_sources`.
996 output_path : `str`, optional
997 Filename to write the combined figure to as a PNG. Parent
998 directories are created if missing.
999 display_figure : `bool`, optional
1000 If True, display the figure inline (notebook).
1001 column_labels : pair of `str`, optional
1002 Labels appended to each cutout's title to identify the column.
1003 figsize_per_row : `float`, optional
1004 Height in inches allocated to each cutout row.
1005 All other kwargs:
1006 Forwarded to the cutout renderer; same meaning as in
1007 `plot_cutouts_with_object_markers`.
1009 Returns
1010 -------
1011 sources, related_objects1, related_objects2 : `pandas.DataFrame`
1012 The catalogs returned by `find_objects_sharing_sources`.
1013 """
1014 import matplotlib.pyplot as plt
1016 if image_type not in image_datasets:
1017 raise ValueError(
1018 f"image_type must be one of {sorted(image_datasets)}, "
1019 f"got {image_type!r}")
1021 sources, ro1, ro2 = find_objects_sharing_sources(
1022 diaObjectId, sources1, sources2, objects1, objects2,
1023 max_distance_arcsec=max_distance_arcsec)
1025 if len(sources) == 0:
1026 print(f"No diaSources in the cluster for "
1027 f"diaObjectId={diaObjectId}")
1028 return sources, ro1, ro2
1030 left_overlays = _prepare_object_overlays(ro1, palette)
1031 right_overlays = _prepare_object_overlays(ro2, palette)
1033 # diaSourceId -> run-2 diaObjectId, so the right panel can color
1034 # each diaSource by its run-2 owner. The left panel uses the
1035 # default (sources["diaObjectId"], the run-1 owner) since
1036 # `sources` is a slice of `sources1`.
1037 src_to_obj2 = dict(zip(
1038 sources2["diaSourceId"].to_numpy(),
1039 sources2["diaObjectId"].to_numpy()))
1040 right_match_ids = [src_to_obj2.get(int(sid))
1041 for sid in sources["diaSourceId"]]
1043 n_rows = len(sources)
1044 # One subfigure per row so each row can carry a single shared
1045 # suptitle above both panels; the per-axes title is then just the
1046 # column label.
1047 fig = plt.figure(figsize=(8, figsize_per_row * n_rows),
1048 constrained_layout=True)
1049 subfigs = np.atleast_1d(fig.subfigures(n_rows, 1, squeeze=False).ravel())
1051 common_kw = dict(
1052 marker_size=marker_size, marker_symbol=marker_symbol,
1053 source_marker_size=source_marker_size,
1054 current_source_marker_size=current_source_marker_size,
1055 current_source_color=current_source_color)
1057 for i, row in enumerate(sources.itertuples(index=False)):
1058 sf = subfigs[i]
1059 sf.suptitle(
1060 f"diaSourceId={int(row.diaSourceId)} "
1061 f"({image_type}, visit={int(row.visit)}, "
1062 f"det={int(row.detector)})",
1063 fontsize="small")
1064 ax_left, ax_right = sf.subplots(1, 2)
1065 # Load once; both panels in this row use the same image.
1066 cutout_data = _load_cutout(butler, row, size=size,
1067 image_type=image_type,
1068 image_datasets=image_datasets)
1069 _render_cutout_axes(ax_left, row, cutout_data, sources,
1070 *left_overlays,
1071 title="", subtitle=column_labels[0],
1072 **common_kw)
1073 _render_cutout_axes(ax_right, row, cutout_data, sources,
1074 *right_overlays,
1075 title="", subtitle=column_labels[1],
1076 source_match_ids=right_match_ids,
1077 **common_kw)
1079 if output_path is not None:
1080 out_dir = os.path.dirname(output_path)
1081 if out_dir:
1082 os.makedirs(out_dir, exist_ok=True)
1083 fig.savefig(output_path, bbox_inches="tight")
1084 if display_figure:
1085 display(fig)
1086 plt.close(fig)
1088 return sources, ro1, ro2
1091def get_xy_from_source_table(table, wcs, degrees=None):
1092 """Convert ra/dec coordinates in an astropy table/pandas data frame to
1093 pixel x/y positions.
1095 Parameters
1096 ----------
1097 table : `astropy.table.Table` or `pandas.DataFrame`
1098 Sources with ``ra``/``dec`` (degrees) or ``coord_ra``/``coord_dec``
1099 (radians) columns.
1100 wcs : `lsst.images.SkyProjection` or `lsst.afw.geom.SkyWcs`
1101 Transform from sky to pixel coordinates.
1102 degrees : `bool`, optional
1103 Whether the coordinates are in degrees; inferred from the column
1104 names if not given.
1106 Returns
1107 -------
1108 xy : `astropy.table.Table`
1109 Pixel positions, in columns ``x`` and ``y``.
1110 """
1111 try:
1112 ra = table['ra']
1113 dec = table['dec']
1114 inferred_degrees = True
1115 except KeyError:
1116 ra = table['coord_ra']
1117 dec = table['coord_dec']
1118 inferred_degrees = False
1119 if degrees is None:
1120 degrees = inferred_degrees
1122 if not degrees:
1123 ra, dec = np.degrees(ra), np.degrees(dec)
1124 x, y = _sky_to_pixel(wcs, ra, dec)
1125 return astropy.table.Table.from_pandas(pd.DataFrame({'x': x, 'y': y}))
1128# Palette used by the `color_by` flag-bucketing mode. Sources with none of
1129# the requested flags set get the residual color "white", which is kept out
1130# of this palette so it never collides with a flagged bucket.
1131_FLAG_PALETTE = ("red", "orange", "yellow", "magenta", "cyan", "green")
1134def _group_sources_by_flag(table, flag_names, palette=_FLAG_PALETTE):
1135 """Split a source table into per-flag buckets for color-coded overlay.
1137 Each row is assigned to the first flag in ``flag_names`` whose column
1138 is True; remaining rows go into a residual "no flag" bucket. Names that
1139 aren't present as columns in ``table`` are silently skipped.
1141 Parameters
1142 ----------
1143 table : table-like
1144 Anything that supports ``len(table)``, ``table[name]`` returning a
1145 boolean-coercible column, and ``table[bool_array]`` row selection.
1146 flag_names : sequence of str
1147 Column names to group on. Order determines color *and* priority
1148 when a row has multiple flags set.
1149 palette : sequence of str, optional
1150 Cycle of display ``ctype`` values to assign in order.
1152 Returns
1153 -------
1154 buckets : list of ``(subset_table, ctype, legend)`` tuples.
1155 """
1156 n = len(table)
1157 if n == 0:
1158 return []
1159 remaining = np.ones(n, dtype=bool)
1160 buckets = []
1161 for i, flag in enumerate(flag_names):
1162 try:
1163 col = table[flag]
1164 except KeyError:
1165 continue
1166 mask = np.asarray(col, dtype=bool) & remaining
1167 if mask.any():
1168 buckets.append((table[mask], palette[i % len(palette)], flag))
1169 remaining = remaining & ~mask
1170 if remaining.any():
1171 buckets.append((table[remaining], "white", "no flag"))
1172 return buckets
1175def _line_segments_from(source_table, wcs, *, flag_col, angle_col, ctype,
1176 length_col=None, fixed_length=None,
1177 min_length=None, length_scale=1.0, thickness=1.5):
1178 """Build centered line-segment endpoints from a source catalog.
1180 Sources are selected by ``flag_col`` (rows where the boolean column
1181 is True) when given, then optionally by ``length_col > min_length``
1182 (only meaningful when ``length_col`` is set), then always by
1183 ``sky_source == False`` when that column is present. Endpoints are
1184 the source centroid ± half the (scaled) length along ``angle_col``.
1186 Exactly one of ``length_col`` (per-source measurement, e.g.
1187 ``ext_trailedSources_Naive_length`` in pixels) or ``fixed_length``
1188 (constant, in pixels — used for markers whose real separation is
1189 below display resolution) must be given. Length units are pixels
1190 and angle units are radians in the detector frame, matching the
1191 raw ``ip_diffim`` and ``ext_trailedSources`` measurement outputs
1192 on ``dia_source_unfiltered``. ``length_scale`` multiplies measured
1193 lengths after filtering; it is a no-op when ``fixed_length`` is
1194 used, since a "fixed" length by definition is not magnified.
1196 ``thickness`` is the stroke width forwarded to ``afw_display.line``;
1197 stored as ``size`` in the returned dict.
1199 Returns a dict ``{"x1", "y1", "x2", "y2", "ctype", "size"}`` of
1200 numpy arrays and scalars, or ``None`` if the table is missing/empty
1201 or any referenced column is absent.
1202 """
1203 if (length_col is None) == (fixed_length is None):
1204 raise ValueError("_line_segments_from: pass exactly one of "
1205 "length_col or fixed_length")
1206 if source_table is None or len(source_table) == 0:
1207 return None
1208 try:
1209 angle = np.asarray(source_table[angle_col], dtype=float)
1210 mask = (np.asarray(source_table[flag_col], dtype=bool) if flag_col is not None
1211 else np.ones(len(source_table), dtype=bool))
1212 if length_col is not None:
1213 length = np.asarray(source_table[length_col], dtype=float)
1214 else:
1215 length = np.full(len(source_table), fixed_length, dtype=float)
1216 except KeyError:
1217 return None
1218 if length_col is not None and min_length is not None:
1219 mask = mask & (length > min_length)
1220 try:
1221 mask = mask & ~np.asarray(source_table["sky_source"], dtype=bool)
1222 except KeyError:
1223 pass
1224 if not mask.any():
1225 return None
1226 xy = get_xy_from_source_table(source_table[mask], wcs)
1227 x0 = xy["x"].data
1228 y0 = xy["y"].data
1229 effective_scale = length_scale if fixed_length is None else 1.0
1230 half = length[mask] * effective_scale / 2.0
1231 dx = half * np.cos(angle[mask])
1232 dy = half * np.sin(angle[mask])
1233 return {"x1": x0 - dx, "y1": y0 - dy,
1234 "x2": x0 + dx, "y2": y0 + dy,
1235 "ctype": ctype, "size": thickness}
1238@dataclass
1239class _OverlayData:
1240 """Everything `display_images` / `display_images_ab` draw on one frame.
1242 ``unfiltered_footprints`` is populated only when the caller asked for
1243 footprint-style rendering of the unfiltered catalog (see
1244 ``unfiltered_as_footprints`` in `_collect_overlays`); it holds
1245 ``(catalog, color, label)`` buckets ready for `_overlay_footprint_layers`,
1246 and in that case the unfiltered ``+`` marker is omitted from
1247 ``overlays``. Otherwise it is ``None`` and the unfiltered catalog is a
1248 normal marker entry in ``overlays``.
1249 """
1250 overlays: list = field(default_factory=list)
1251 reliability_labels: dict | None = None
1252 solar_system_labels: dict | None = None
1253 dipole_segments: dict | None = None
1254 trail_segments: dict | None = None
1255 unfiltered_footprints: list | None = None
1258def _collect_overlays(butler, data_id, wcs, *,
1259 reliability_threshold,
1260 show_unfiltered, show_trailed,
1261 show_rejected, show_standardized, show_marginal,
1262 show_kernel_sources, show_solar_system, show_apdb,
1263 show_reliability_labels, show_dipoles,
1264 show_trail_geometry, line_length_scale,
1265 color_by, unfiltered_as_footprints=False):
1266 """Load catalogs from one butler and build the overlay record list.
1268 Shared between `display_images` and `display_images_ab`. Catalogs that
1269 aren't present for this dataId are silently skipped.
1271 When ``unfiltered_as_footprints`` is True the ``dia_source_unfiltered``
1272 catalog is *not* added to ``overlays`` as ``+`` markers; instead its
1273 non-sky rows are returned as `_OverlayData.unfiltered_footprints`
1274 ``(catalog, color, label)`` buckets for `_overlay_footprint_layers` to
1275 draw as Firefly footprint layers. With ``color_by`` the buckets are the
1276 flag partition (`_group_sources_by_flag`); otherwise a single red
1277 bucket. The dipole/trail line segments still come from the same
1278 unfiltered catalog either way.
1280 Returns
1281 -------
1282 data : `_OverlayData`
1283 ``overlays`` is a list of ``(x_arr, y_arr, symbol, size, ctype,
1284 legend)`` marker tuples; ``reliability_labels`` /
1285 ``solar_system_labels`` are ``{"x", "y", ...}`` dicts (or None) for
1286 text annotations; ``dipole_segments`` / ``trail_segments`` are
1287 line-segment dicts (or None); ``unfiltered_footprints`` is the
1288 footprint-bucket list described above (or None).
1289 """
1290 def _try_get(dataset):
1291 try:
1292 return butler.get(dataset, data_id)
1293 except DatasetNotFoundError:
1294 return None
1296 overlays = []
1298 def _add(table, *, symbol, size, ctype, legend, use_radec=True):
1299 if table is None or len(table) == 0:
1300 return
1301 if use_radec:
1302 xy = get_xy_from_source_table(table, wcs)
1303 x_arr = xy["x"].data
1304 y_arr = xy["y"].data
1305 else:
1306 x_arr = table["x"].data
1307 y_arr = table["y"].data
1308 overlays.append((x_arr, y_arr, symbol, size, ctype, legend))
1310 # Catalogs are added in AP-pipeline-creation order
1311 if show_kernel_sources:
1312 _add(_try_get("difference_kernel_sources"),
1313 symbol="o", size=12, ctype="green",
1314 legend="psf-matching kernel source")
1316 # Load dia_source_unfiltered once — it backs the unfiltered marker
1317 # overlay AND the dipole/trail line-segment overlays below (which
1318 # always read from the unfiltered catalog because dipoles and long
1319 # trails are filtered out of the downstream catalogs).
1320 unfiltered = None
1321 if show_unfiltered or show_dipoles or show_trail_geometry:
1322 unfiltered = _try_get("dia_source_unfiltered")
1324 unfiltered_footprints = None
1325 if show_unfiltered and unfiltered is not None and len(unfiltered) > 0:
1326 non_sky = unfiltered[~unfiltered["sky_source"]]
1327 if unfiltered_as_footprints:
1328 # Draw as footprints instead of markers: build one
1329 # (catalog, color, label) bucket per color for
1330 # `_overlay_footprint_layers`. color_by splits by flag
1331 # (deterministic); otherwise a single red bucket.
1332 if color_by:
1333 unfiltered_footprints = [
1334 (sub, ctype, f"unfiltered: {flag}")
1335 for sub, ctype, flag
1336 in _group_sources_by_flag(non_sky, color_by)]
1337 else:
1338 unfiltered_footprints = [(non_sky, "red", "unfiltered candidate")]
1339 elif color_by:
1340 for sub, ctype, flag in _group_sources_by_flag(non_sky, color_by):
1341 _add(sub, symbol="+", size=10, ctype=ctype,
1342 legend=f"unfiltered: {flag}")
1343 else:
1344 _add(non_sky, symbol="+", size=10, ctype="red",
1345 legend="unfiltered candidate")
1346 if show_rejected:
1347 _add(_try_get("rejected_dia_source"),
1348 symbol="+", size=10, ctype="orange", legend="rejected diaSource")
1349 if show_trailed:
1350 _add(_try_get("long_trailed_source_detector"),
1351 symbol="x", size=30, ctype="magenta", legend="long-trailed source")
1353 # Stash the standardized catalog + projected xy for reuse by the
1354 # geometry overlays below. `standardizeDiaSource` runs between
1355 # filterDiaSource and associateApdb; when the pipeline stops before
1356 # the APDB ingest, dia_source_detector is the last diaSource catalog
1357 # available.
1358 standardized_data = None
1359 if show_standardized:
1360 standardized = _try_get("dia_source_detector")
1361 if standardized is not None and len(standardized) > 0:
1362 xy = get_xy_from_source_table(standardized, wcs)
1363 x_arr = xy["x"].data
1364 y_arr = xy["y"].data
1365 overlays.append((x_arr, y_arr, "+", 10, "blue", "standardized diaSource"))
1366 standardized_data = {"catalog": standardized, "x": x_arr, "y": y_arr}
1368 # Load dia_source_apdb once: it backs the APDB reliability overlay and
1369 # also supplies pixel x/y for the solar-system overlay (ss_source_detector
1370 # carries only the matched diaSourceId, not coordinates).
1371 dia_apdb = None
1372 if show_solar_system or show_apdb:
1373 dia_apdb = _try_get("dia_source_apdb")
1375 if show_apdb and dia_apdb is not None and len(dia_apdb) > 0:
1376 good_mask = dia_apdb["reliability"] > reliability_threshold
1377 _add(dia_apdb[good_mask], symbol="o", size=14, ctype="blue", use_radec=False,
1378 legend=f"APDB, reliability > {reliability_threshold:g}")
1379 _add(dia_apdb[~good_mask], symbol="o", size=14, ctype="red", use_radec=False,
1380 legend=f"APDB, reliability <= {reliability_threshold:g}")
1382 solar_system_labels = None
1383 if show_solar_system:
1384 ss = _try_get("ss_source_detector")
1385 if (ss is not None and len(ss) > 0
1386 and dia_apdb is not None and len(dia_apdb) > 0):
1387 # ss_source_detector lacks coords; match each diaSourceId to
1388 # the APDB row to recover its pixel x/y.
1389 ss_ids = np.asarray(ss["diaSourceId"])
1390 apdb_ids = np.asarray(dia_apdb["diaSourceId"])
1391 idx_in_apdb = pd.Series(np.arange(len(apdb_ids)), index=apdb_ids).reindex(ss_ids)
1392 keep = idx_in_apdb.notna().to_numpy()
1393 if keep.any():
1394 apdb_idx = idx_in_apdb.dropna().astype(int).to_numpy()
1395 x_arr = np.asarray(dia_apdb["x"])[apdb_idx]
1396 y_arr = np.asarray(dia_apdb["y"])[apdb_idx]
1397 designation = np.asarray(ss["designation"])[keep]
1398 overlays.append((x_arr, y_arr, "o", 16, "cyan", "solar-system match"))
1399 solar_system_labels = {"x": x_arr, "y": y_arr, "designation": designation, }
1401 if show_marginal:
1402 _add(_try_get("marginal_new_dia_source"),
1403 symbol="+", size=10, ctype="yellow", legend="marginal new diaSource")
1405 # Reliability text is drawn at most once per diaSource: prefer APDB
1406 # good sources when APDB is displayed, and otherwise annotate every
1407 # standardized row (they've been pipeline-filtered to high
1408 # reliability already). The unfiltered catalog doesn't carry a final
1409 # reliability score, so it never provides labels.
1410 reliability_labels = None
1411 apdb_shown = show_apdb and dia_apdb is not None and len(dia_apdb) > 0
1412 if show_reliability_labels:
1413 if apdb_shown:
1414 rel = np.asarray(dia_apdb["reliability"])
1415 mask = rel > reliability_threshold
1416 if mask.any():
1417 reliability_labels = {"x": np.asarray(dia_apdb["x"])[mask],
1418 "y": np.asarray(dia_apdb["y"])[mask],
1419 "reliability": rel[mask]}
1420 elif standardized_data is not None:
1421 reliability_labels = {
1422 "x": standardized_data["x"],
1423 "y": standardized_data["y"],
1424 "reliability": np.asarray(standardized_data["catalog"]["reliability"]),
1425 }
1427 # Dipole and trail line segments both come from dia_source_unfiltered:
1428 # it's the earliest AP-pipeline catalog and thus a superset of the
1429 # downstream ones that get dipoles/long trails filtered out, and it
1430 # carries the raw ip_diffim / ext_trailedSources measurement columns
1431 # in native pixel + detector-radian units — exactly the coordinate
1432 # system the endpoint math lives in.
1433 dipole_segments = None
1434 if show_dipoles:
1435 # Fixed 10 px length: the measured ``ip_diffim_DipoleFit_separation``
1436 # is sub-pixel for the vast majority of classified dipoles (median
1437 # ~0.09 px in typical data), so drawing at the measured length
1438 # would hide them. The line here is a fixed-size *marker* of the
1439 # dipole's orientation, not a physical extent.
1440 dipole_segments = _line_segments_from(
1441 unfiltered, wcs,
1442 flag_col="ip_diffim_DipoleFit_classification",
1443 angle_col="ip_diffim_DipoleFit_orientation",
1444 ctype="white", fixed_length=10.0, thickness=3.0)
1446 # 3 px threshold: below this the trail measurement is dominated by
1447 # noise on point-like sources. Long-trailed sources removed by
1448 # filterDiaSource still show up via ``show_trailed`` as ``x`` markers.
1449 trail_segments = None
1450 if show_trail_geometry:
1451 trail_segments = _line_segments_from(
1452 unfiltered, wcs, flag_col=None,
1453 length_col="ext_trailedSources_Naive_length",
1454 angle_col="ext_trailedSources_Naive_angle",
1455 ctype="magenta", min_length=3.0,
1456 length_scale=line_length_scale)
1458 return _OverlayData(
1459 overlays=overlays,
1460 reliability_labels=reliability_labels,
1461 solar_system_labels=solar_system_labels,
1462 dipole_segments=dipole_segments,
1463 trail_segments=trail_segments,
1464 unfiltered_footprints=unfiltered_footprints,
1465 )
1468def _print_overlay_legend(overlays, header, indent=""):
1469 """Print a one-line-per-overlay legend for a single panel."""
1470 print(f"{indent}{header}")
1471 for x_arr, _, symbol, _, ctype, legend in overlays:
1472 print(f"{indent} {len(x_arr):5d} {ctype:>8s} {symbol} {legend}")
1475def _print_footprint_legend(buckets, indent=""):
1476 """Print a one-line-per-color summary of the unfiltered footprint layers.
1478 ``buckets`` is the `_OverlayData.unfiltered_footprints` list of
1479 ``(catalog, color, label)`` triples; does nothing when it is empty or
1480 None (i.e. the unfiltered catalog was drawn as markers).
1481 """
1482 if not buckets:
1483 return
1484 total = sum(len(cat) for cat, _, _ in buckets)
1485 print(f"{indent}{total:5d} footprints unfiltered ({len(buckets)} colors):")
1486 for cat, color, label in buckets:
1487 print(f"{indent} {len(cat):5d} {color:>11s} {label}")
1490def _resolve_unfiltered_footprints(unfiltered_style, backend, color_by):
1491 """Decide whether to draw the unfiltered catalog as footprints, warning
1492 about the caveats.
1494 Footprints need Firefly's native overlay, so a non-firefly backend
1495 falls back to ``+`` markers. When footprints are combined with
1496 ``color_by``, stale layers from a previous call are not auto-erased
1497 (Firefly exposes no API to delete or enumerate footprint layers), so
1498 warn that clearing them before re-running is the caller's
1499 responsibility.
1500 """
1501 if unfiltered_style not in ("footprint", "marker"):
1502 raise ValueError("unfiltered_style must be 'footprint' or 'marker', "
1503 f"got {unfiltered_style!r}")
1504 use_footprints = unfiltered_style == "footprint" and backend == "firefly"
1505 if unfiltered_style == "footprint" and backend != "firefly":
1506 print(f"WARNING: unfiltered_style='footprint' needs the 'firefly' "
1507 f"backend; falling back to '+' markers for backend={backend!r}.")
1508 if use_footprints and color_by:
1509 print("WARNING: color_by footprint layers are not auto-erased; "
1510 "re-running may leave stale footprints from a previous call. "
1511 "Clear the frame's footprint layers before re-running.")
1512 return use_footprints
1515def _draw_line_segments(afw_display, segments):
1516 """Draw one batch of centered line segments on the active frame."""
1517 if segments is None:
1518 return
1519 ctype = segments["ctype"]
1520 size = segments["size"]
1521 for x1, y1, x2, y2 in zip(segments["x1"], segments["y1"],
1522 segments["x2"], segments["y2"]):
1523 afw_display.line([(float(x1), float(y1)), (float(x2), float(y2))],
1524 ctype=ctype, size=size)
1527def _draw_overlays_on_current_frame(afw_display, overlays,
1528 reliability_labels, solar_system_labels,
1529 dipole_segments=None,
1530 trail_segments=None,
1531 label_size=3):
1532 """Stamp one set of overlays + optional reliability and solar-system
1533 designation labels onto the active frame.
1535 ``label_size`` is the text size (in pixels) used for both label sets.
1536 ``dipole_segments`` and ``trail_segments`` are optional line-segment
1537 dicts (see `_line_segments_from`).
1538 """
1539 # Scale the text offset with the size so larger labels still clear the
1540 # circle markers they annotate.
1541 label_offset = max(14, 2 * label_size)
1542 with afw_display.Buffering():
1543 for x_arr, y_arr, symbol, size, ctype, _ in overlays:
1544 for x, y in zip(x_arr, y_arr):
1545 afw_display.dot(symbol, x, y, size=size, ctype=ctype)
1546 # Trails first, dipoles on top: a source flagged as a dipole is
1547 # the more actionable pipeline-quality issue, so it wins any
1548 # pixel overlap with the trail line.
1549 _draw_line_segments(afw_display, trail_segments)
1550 _draw_line_segments(afw_display, dipole_segments)
1551 if reliability_labels is not None:
1552 # Offset the score text so it doesn't sit on top of the marker.
1553 for r, x, y in zip(reliability_labels["reliability"],
1554 reliability_labels["x"],
1555 reliability_labels["y"]):
1556 afw_display.dot(f"{r:.2f}", x + label_offset, y,
1557 size=label_size, ctype="cyan")
1558 if solar_system_labels is not None:
1559 # Offset SS designations *below* the marker so they don't
1560 # overplot any reliability score drawn to the right.
1561 for desig, x, y in zip(solar_system_labels["designation"],
1562 solar_system_labels["x"],
1563 solar_system_labels["y"]):
1564 afw_display.dot(str(desig), x, y + label_offset,
1565 size=label_size, ctype="cyan")
1568def _get_image(butler, image_datasets, image_type, data_id):
1569 """Read one AP image in `lsst.images` form.
1571 Parameters
1572 ----------
1573 butler : `lsst.daf.butler.Butler`
1574 Butler to read from.
1575 image_datasets : `dict` [`str`, `str`]
1576 Mapping from image type to dataset name.
1577 image_type : {"science", "template", "difference"}
1578 Image to read.
1579 data_id : `dict`
1580 Data ID of the image.
1582 Returns
1583 -------
1584 image : `lsst.images.VisitImage`
1585 The image; a legacy template is given the science image's
1586 detector, which it lacks.
1587 """
1588 detector = None
1589 if image_type == "template": 1589 ↛ 1590line 1589 didn't jump to line 1590 because the condition on line 1589 was never true
1590 detector = butler.get(f"{image_datasets['science']}.detector", data_id)
1591 return get_visit_image(butler, image_datasets[image_type], data_id, image_type=image_type,
1592 detector=detector)
1595def _science_bbox(butler, image_datasets, data_id):
1596 """Read the bounding box of a science image as an `lsst.images.Box`."""
1597 bbox = butler.get(f"{image_datasets['science']}.bbox", data_id)
1598 return Box.from_legacy(bbox) if isinstance(bbox, lsst.geom.Box2I) else bbox
1601def _legacy_for_display(image):
1602 """Return the legacy form of an image, which afw displays take."""
1603 return image.to_legacy() if hasattr(image, "to_legacy") else image
1606def _show_image(afw_display, image, title, strip_metadata):
1607 """Show an image on the current frame of an afw display.
1609 Parameters
1610 ----------
1611 afw_display : `lsst.afw.display.Display`
1612 Display to draw on.
1613 image : `lsst.images.MaskedImage`, `lsst.images.cells.CellCoadd`, or \
1614 `lsst.afw.image.Exposure`
1615 Image to show; afw displays take only the legacy form, so an
1616 `lsst.images` image is converted.
1617 title : `str`
1618 Frame title.
1619 strip_metadata : `bool`
1620 Drop ``LTV1``/``LTV2`` from the shown exposure's metadata.
1622 Returns
1623 -------
1624 exposure : `lsst.afw.image.Exposure`
1625 The exposure sent to the display.
1626 """
1627 exposure = _legacy_for_display(image)
1628 if strip_metadata:
1629 _strip_ds9_metadata(exposure)
1630 afw_display.image(exposure, title=title)
1631 return exposure
1634def _strip_ds9_metadata(*exposures):
1635 """Drop LTV1/LTV2 keys from each exposure's metadata in place."""
1636 for exp in exposures:
1637 md = exp.metadata
1638 for k in ("LTV1", "LTV2"):
1639 if md.exists(k):
1640 md.remove(k)
1643def _erase_current_frame_regions(afw_display):
1644 """Clear all region markers on the currently-selected frame.
1646 The firefly backend caches ``_regionLayerId`` on its impl and only
1647 refreshes it inside ``_flush()``. Calling ``erase()`` right after
1648 switching frames therefore issues ``delete_region_layer`` with the
1649 *previous* frame's layer id under the current frame's plot id, so
1650 nothing gets removed. Refreshing the cached id from the current
1651 frame before erasing fixes it. Backends that don't cache a layer
1652 id (e.g. ds9) fall through to a plain ``erase()``.
1653 """
1654 impl = getattr(afw_display, "_impl", None)
1655 if impl is not None and hasattr(impl, "_getRegionLayerId"):
1656 impl._regionLayerId = impl._getRegionLayerId()
1657 afw_display.erase()
1660def display_images(butler, visit, detector, backend="firefly", *,
1661 reliability_threshold=0.1,
1662 show_unfiltered=True,
1663 show_trailed=True,
1664 show_rejected=True,
1665 show_standardized=True,
1666 show_marginal=True,
1667 show_kernel_sources=True,
1668 show_solar_system=True,
1669 show_apdb=True,
1670 show_reliability_labels=True,
1671 show_dipoles=True,
1672 show_trail_geometry=True,
1673 line_length_scale=1.0,
1674 label_size=3,
1675 color_by=None,
1676 unfiltered_style="footprint",
1677 mask_transparency=80,
1678 strip_metadata=True,
1679 skymap=None,
1680 skymap_ctype="green",
1681 skymap_label_size=1.5,
1682 image_datasets=_IMAGE_DATASETS,
1683 use_fakes=False,
1684 dry_run=False):
1685 """Display the science, template, and difference images for a given
1686 visit+detector with diagnostic catalog markers overlaid.
1688 Three frames are produced (science, template, difference) and the same
1689 overlays are drawn on each. Catalogs that are missing from the butler
1690 are silently skipped, so the same call works against partial outputs.
1692 Default overlay key. Rows are in AP-pipeline creation order, so the
1693 last marker drawn at any pixel reflects the latest classification the
1694 pipeline assigned. Circle sizes step by 2 so successive ``o`` markers
1695 nest rather than stack.
1697 ============================= ========= ==== ==========================
1698 catalog symbol size color
1699 ============================= ========= ==== ==========================
1700 psf-matching kernel sources ``o`` 12 green
1701 unfiltered candidates footprint -- red (``unfiltered_style``)
1702 rejected diaSources ``+`` 10 orange
1703 long-trailed sources ``x`` 30 magenta
1704 standardized diaSources ``+`` 10 blue
1705 APDB, reliability > threshold ``o`` 14 blue (+ score text)
1706 APDB, reliability ≤ threshold ``o`` 14 red
1707 solar-system matches ``o`` 16 cyan
1708 marginal new diaSources ``+`` 10 yellow
1709 ============================= ========= ==== ==========================
1711 By default the ``dia_source_unfiltered`` catalog is drawn as Firefly
1712 footprint outlines rather than ``+`` markers; see ``unfiltered_style``.
1714 Parameters
1715 ----------
1716 butler : `lsst.daf.butler.Butler`
1717 Butler to load data from.
1718 visit, detector : `int`
1719 Visit and detector ids to load data for.
1720 backend : `str`, optional
1721 afw display backend (typically "firefly" or "ds9").
1722 reliability_threshold : `float`, optional
1723 APDB diaSources with reliability strictly greater than this are
1724 drawn as "good" (blue); the rest as "bad" (red).
1725 show_unfiltered, show_trailed, show_rejected, show_standardized,
1726 show_marginal, show_kernel_sources, show_solar_system,
1727 show_apdb : `bool`, optional
1728 Toggle individual catalog overlays. ``show_kernel_sources``
1729 loads ``difference_kernel_sources``, the PSF-matching constraint
1730 sources from image subtraction — useful for seeing where the
1731 kernel was actually anchored vs extrapolated.
1732 show_reliability_labels : `bool`, optional
1733 If True, annotate each good APDB diaSource with its reliability score.
1734 show_dipoles : `bool`, optional
1735 If True, draw a 10-px white line segment through each source in
1736 ``dia_source_unfiltered`` with ``ip_diffim_DipoleFit_classification``
1737 set, oriented along ``ip_diffim_DipoleFit_orientation`` (radians,
1738 detector-frame). The length is fixed rather than measured because
1739 ``ip_diffim_DipoleFit_separation`` is sub-pixel for the vast
1740 majority of classified dipoles; the line is a fixed-size marker
1741 of orientation rather than a physical extent. Sourced from the
1742 unfiltered catalog rather than the standardized or APDB catalogs
1743 because ``filterDiaSource`` removes dipoles before those stages.
1744 show_trail_geometry : `bool`, optional
1745 If True, draw a magenta line segment along
1746 ``ext_trailedSources_Naive_angle`` with length
1747 ``ext_trailedSources_Naive_length`` for every source in
1748 ``dia_source_unfiltered`` whose trail length exceeds 3 px.
1749 Long-trailed sources removed by ``filterDiaSource`` still show
1750 up separately under ``show_trailed`` as ``x`` markers.
1751 line_length_scale : `float`, optional
1752 Multiplicative factor applied to the drawn length of the trail
1753 line segments *after* the 3 px trail filter, so a below-threshold
1754 trail stays hidden regardless of the scale. Does not affect the
1755 dipole marker, which is drawn at a fixed 10 px length by design.
1756 Default 1.0 (draw trails at their measured length); use larger
1757 values to make short trails easier to see against the image.
1758 label_size : `int`, optional
1759 Text size (in pixels) for the reliability score and solar-system
1760 designation annotations.
1761 color_by : sequence of `str`, optional
1762 Flag column names from ``dia_source_unfiltered``. When supplied,
1763 the unfiltered-candidate overlay is split into buckets colored by
1764 which named flag fires first (list order = color *and* priority),
1765 with a residual white bucket for rows that match none of them.
1766 Unknown column names are silently skipped. Applies whether the
1767 unfiltered catalog is drawn as footprints or markers. Example::
1769 color_by=["pixelFlags_bad", "pixelFlags_edge",
1770 "ip_diffim_DipoleFit_classification",
1771 "pixelFlags_saturated"]
1772 unfiltered_style : {"footprint", "marker"}, optional
1773 How to draw ``dia_source_unfiltered``. ``"footprint"`` (default)
1774 overlays each source's `Footprint` outline using Firefly's native
1775 footprint rendering; ``"marker"`` draws the old red ``+`` markers.
1776 Footprints need the ``"firefly"`` backend, so any other backend
1777 (e.g. ds9) silently falls back to markers. Two caveats with
1778 footprints, both warned about at call time: (1) when combined with
1779 ``color_by``, footprint layers from a previous call are *not*
1780 auto-erased (Firefly offers no way to delete or enumerate them), so
1781 clearing them before re-running is your responsibility; (2)
1782 switching ``unfiltered_style`` from ``"footprint"`` to ``"marker"``
1783 leaves the previous footprints on the frame — re-run in footprint
1784 mode or reset the frame to clear them. Re-running in the default
1785 (no ``color_by``) footprint mode overwrites the single layer in
1786 place, so it is unaffected.
1787 mask_transparency : `int` or `None`, optional
1788 Mask-plane transparency forwarded to the display (0 = opaque,
1789 100 = fully transparent). Pass ``None`` to leave the backend's
1790 current setting untouched.
1791 strip_metadata : `bool`, optional
1792 Drop ``LTV1``/``LTV2`` keywords from each exposure's metadata
1793 before sending to the backend. Needed for ds9 to align frames.
1794 skymap : `lsst.skymap.BaseSkyMap`, optional
1795 If supplied, overlay the boundaries of every tract/patch that
1796 touches each frame, labeled ``tract,patch``.
1797 skymap_ctype : `str`, optional
1798 Display color for the tract/patch outlines and labels.
1799 skymap_label_size : `float`, optional
1800 Text size for the ``tract,patch`` labels.
1801 image_datasets : `dict` [`str`, `str`], optional
1802 Mapping from image-type key (``"science"``, ``"template"``,
1803 ``"difference"``) to butler dataset name. Override to point at
1804 alternate dataset types.
1805 use_fakes : `bool`, optional
1806 If True, load the fake-source-injected versions of the science
1807 and template images: ``fakes_`` prefix on the science dataset
1808 and ``injectedTemplate_`` prefix on the template dataset. The
1809 difference image and every catalog keep their non-prefixed
1810 names, per the fake-source pipeline's output convention.
1811 Default False.
1812 dry_run : `bool`, optional
1813 If True, load every requested image and catalog and print the
1814 overlay/footprint legends, but skip constructing the afw display
1815 and drawing anything. Useful for sanity-checking which datasets
1816 are available for a (visit, detector) without opening a viewer.
1817 Default False.
1818 """
1819 data_id = {"visit": visit, "detector": detector}
1820 image_datasets = _apply_fakes_prefix(image_datasets, use_fakes)
1821 use_footprints = _resolve_unfiltered_footprints(
1822 unfiltered_style, backend, color_by)
1824 diffim = _get_image(butler, image_datasets, "difference", data_id)
1825 science = _get_image(butler, image_datasets, "science", data_id)
1826 template = get_visit_image(butler, image_datasets["template"], data_id, image_type="template",
1827 detector=science.detector)
1828 template = template[science.bbox]
1829 images = {"science": science, "template": template, "difference": diffim}
1831 data = _collect_overlays(
1832 butler, data_id, diffim.sky_projection,
1833 reliability_threshold=reliability_threshold,
1834 show_unfiltered=show_unfiltered,
1835 show_trailed=show_trailed, show_rejected=show_rejected,
1836 show_standardized=show_standardized,
1837 show_marginal=show_marginal,
1838 show_kernel_sources=show_kernel_sources,
1839 show_solar_system=show_solar_system,
1840 show_apdb=show_apdb,
1841 show_reliability_labels=show_reliability_labels,
1842 show_dipoles=show_dipoles,
1843 show_trail_geometry=show_trail_geometry,
1844 line_length_scale=line_length_scale,
1845 color_by=color_by,
1846 unfiltered_as_footprints=use_footprints,
1847 )
1848 _print_overlay_legend(
1849 data.overlays, f"visit={visit}, detector={detector} -- overlay legend:")
1850 _print_footprint_legend(data.unfiltered_footprints, indent=" ")
1852 if dry_run:
1853 return
1855 afw_display = lsst.afw.display.Display(backend=backend)
1856 if mask_transparency is not None:
1857 afw_display.setMaskTransparency(mask_transparency)
1858 for frame, image_name in enumerate(("science", "template", "difference")):
1859 afw_display.frame = frame
1860 # Wipe any markers left over from a previous call — `image()`
1861 # only replaces the pixel data, region overlays persist otherwise.
1862 _erase_current_frame_regions(afw_display)
1863 exposure = _show_image(afw_display, images[image_name], image_name, strip_metadata)
1864 _draw_overlays_on_current_frame(
1865 afw_display, data.overlays, data.reliability_labels,
1866 data.solar_system_labels,
1867 dipole_segments=data.dipole_segments,
1868 trail_segments=data.trail_segments,
1869 label_size=label_size)
1870 if data.unfiltered_footprints:
1871 # Per-frame layer prefix so the same catalog drawn on all three
1872 # frames gets distinct Firefly layer ids.
1873 _overlay_footprint_layers(
1874 afw_display, data.unfiltered_footprints,
1875 style="outline", layer_prefix=f"{image_name} unfiltered")
1876 if skymap is not None:
1877 draw_skymap_outlines_afw(afw_display, skymap, exposure.wcs, exposure.getBBox(),
1878 ctype=skymap_ctype, label_size=skymap_label_size)
1880 try:
1881 afw_display.alignImages(match_type="Pixel")
1882 except NotImplementedError:
1883 print(f"WARNING: cannot automatically align and lock images with backend={backend!r}.")
1886def display_images_ab(butler_a, butler_b, visit, detector, *,
1887 image_type="difference",
1888 labels=("A", "B"),
1889 backend="firefly",
1890 reliability_threshold=0.1,
1891 show_unfiltered=True,
1892 show_trailed=True,
1893 show_rejected=True,
1894 show_standardized=True,
1895 show_marginal=True,
1896 show_kernel_sources=True,
1897 show_solar_system=True,
1898 show_apdb=True,
1899 show_reliability_labels=True,
1900 show_dipoles=True,
1901 show_trail_geometry=True,
1902 line_length_scale=1.0,
1903 label_size=3,
1904 color_by=None,
1905 unfiltered_style="footprint",
1906 mask_transparency=80,
1907 strip_metadata=True,
1908 skymap=None,
1909 skymap_ctype="green",
1910 skymap_label_size=1.5,
1911 image_datasets=_IMAGE_DATASETS,
1912 use_fakes=False,
1913 dry_run=False):
1914 """Display one image type side-by-side from two butlers, with overlays.
1916 Loads the same (visit, detector) from ``butler_a`` and ``butler_b``,
1917 places them in two frames, and draws each butler's catalog overlays on
1918 its own frame. Intended for A/B-testing pipeline-config changes that affect
1919 detection or subtraction quality.
1921 Parameters
1922 ----------
1923 butler_a, butler_b : `lsst.daf.butler.Butler`
1924 Two butlers, typically from different pipeline runs of the same data.
1925 visit, detector : `int`
1926 Visit and detector ids to load data for.
1927 image_type : {"science", "template", "difference"}, optional
1928 Which image dataset to compare. Default ``"difference"``.
1929 labels : pair of `str`, optional
1930 Short tags for the two frames; appear in the image title and the
1931 legend header. Default ``("A", "B")``.
1932 backend : `str`, optional
1933 afw display backend (typically "firefly" or "ds9").
1934 reliability_threshold, show_unfiltered, show_trailed, show_rejected,
1935 show_standardized, show_marginal, show_kernel_sources,
1936 show_solar_system, show_apdb, show_reliability_labels, show_dipoles,
1937 show_trail_geometry, line_length_scale, label_size, color_by,
1938 unfiltered_style, mask_transparency, strip_metadata, skymap,
1939 skymap_ctype, skymap_label_size, image_datasets, use_fakes, dry_run
1940 Same meaning as in `display_images`. Applied to both frames; the
1941 tract/patch overlay uses each frame's own exposure WCS. Each frame's
1942 unfiltered footprints get their own per-frame Firefly layers (keyed
1943 by ``labels``), so the two frames don't share layers.
1944 """
1945 if image_type not in image_datasets:
1946 raise ValueError(
1947 f"image_type must be one of {sorted(image_datasets)}, got {image_type!r}")
1948 image_datasets = _apply_fakes_prefix(image_datasets, use_fakes)
1949 data_id = {"visit": visit, "detector": detector}
1951 image_a = _get_image(butler_a, image_datasets, image_type, data_id)
1952 image_b = _get_image(butler_b, image_datasets, image_type, data_id)
1953 if image_type == "template":
1954 # Templates are usually larger than the science footprint; clip them
1955 # to the science bbox so the two frames have matching extents.
1956 image_a = image_a[_science_bbox(butler_a, image_datasets, data_id)]
1957 image_b = image_b[_science_bbox(butler_b, image_datasets, data_id)]
1959 use_footprints = _resolve_unfiltered_footprints(
1960 unfiltered_style, backend, color_by)
1961 common = dict(
1962 reliability_threshold=reliability_threshold,
1963 show_unfiltered=show_unfiltered,
1964 show_trailed=show_trailed, show_rejected=show_rejected,
1965 show_standardized=show_standardized,
1966 show_marginal=show_marginal,
1967 show_kernel_sources=show_kernel_sources,
1968 show_solar_system=show_solar_system,
1969 show_apdb=show_apdb, show_reliability_labels=show_reliability_labels,
1970 show_dipoles=show_dipoles,
1971 show_trail_geometry=show_trail_geometry,
1972 line_length_scale=line_length_scale,
1973 color_by=color_by,
1974 unfiltered_as_footprints=use_footprints,
1975 )
1976 data_a = _collect_overlays(butler_a, data_id, image_a.sky_projection, **common)
1977 data_b = _collect_overlays(butler_b, data_id, image_b.sky_projection, **common)
1979 label_a, label_b = labels
1980 print(f"visit={visit}, detector={detector}: A/B comparison of {image_type!r}")
1981 _print_overlay_legend(data_a.overlays, f"-- {label_a} overlay legend:", indent=" ")
1982 _print_footprint_legend(data_a.unfiltered_footprints, indent=" ")
1983 _print_overlay_legend(data_b.overlays, f"-- {label_b} overlay legend:", indent=" ")
1984 _print_footprint_legend(data_b.unfiltered_footprints, indent=" ")
1986 if dry_run:
1987 return
1989 afw_display = lsst.afw.display.Display(backend=backend)
1990 if mask_transparency is not None:
1991 afw_display.setMaskTransparency(mask_transparency)
1992 for frame, (tag, image, data) in enumerate((
1993 (label_a, image_a, data_a),
1994 (label_b, image_b, data_b))):
1995 afw_display.frame = frame
1996 # Wipe any markers left over from a previous call — `image()`
1997 # only replaces the pixel data, region overlays persist otherwise.
1998 _erase_current_frame_regions(afw_display)
1999 exposure = _show_image(afw_display, image, f"{image_type} ({tag})", strip_metadata)
2000 _draw_overlays_on_current_frame(afw_display, data.overlays,
2001 data.reliability_labels,
2002 data.solar_system_labels,
2003 dipole_segments=data.dipole_segments,
2004 trail_segments=data.trail_segments,
2005 label_size=label_size)
2006 if data.unfiltered_footprints:
2007 # Per-frame layer prefix (the A/B tag) so the two frames'
2008 # footprints get distinct Firefly layer ids.
2009 _overlay_footprint_layers(
2010 afw_display, data.unfiltered_footprints,
2011 style="outline", layer_prefix=f"{tag} unfiltered")
2012 if skymap is not None:
2013 draw_skymap_outlines_afw(afw_display, skymap, exposure.wcs, exposure.getBBox(),
2014 ctype=skymap_ctype, label_size=skymap_label_size)
2016 try:
2017 afw_display.alignImages(match_type="Pixel")
2018 except NotImplementedError:
2019 print(f"WARNING: cannot automatically align and lock images with backend={backend!r}.")
2022@dataclass
2023class _PatchCoverage:
2024 """One coadd patch that contributed to a template."""
2026 tract: int
2027 patch: int
2028 ref: Any
2029 """`lsst.daf.butler.DatasetRef` of the coadd, or None if not found."""
2030 corners_xy: list
2031 """Patch outline projected into difference-image pixels."""
2032 overlap: float
2033 """Fraction of the difference image this patch covers."""
2034 color: str
2035 frame: int | None
2036 """Display frame showing this coadd; None when the coadd is missing."""
2039def _coadd_input_patches(butler, template_dataset, data_id):
2040 """The coadd patches recorded as a template's inputs.
2042 Read as a butler component, so the template's pixels are never
2043 loaded.
2045 Parameters
2046 ----------
2047 butler : `lsst.daf.butler.Butler`
2048 Butler to load from.
2049 template_dataset : `str`
2050 Dataset name of the template (e.g. ``"template_detector"``).
2051 data_id : `dict`
2052 Data id of the template.
2054 Returns
2055 -------
2056 patches : `dict` [`tuple` [`int`, `int`], `lsst.afw.table.ExposureRecord`]
2057 Exposure record of each contributing patch, keyed on
2058 ``(tract, patch)`` and sorted. Empty if the template carries no
2059 coadd inputs.
2060 """
2061 coadd_inputs = butler.get(f"{template_dataset}.coaddInputs", data_id)
2062 ccds = None if coadd_inputs is None else coadd_inputs.ccds
2063 if ccds is None or len(ccds) == 0:
2064 return {}
2065 if not {"tract", "patch"} <= ccds.schema.getNames():
2066 raise RuntimeError(
2067 f"{template_dataset}.coaddInputs has no 'tract'/'patch' fields, so "
2068 "the contributing coadds cannot be identified; this template "
2069 "predates GetTemplateTask recording its input patches.")
2070 patches = {}
2071 for record in ccds:
2072 # Patch ids repeat across tracts, so key on the pair. One record
2073 # per patch is expected; keep the first if that ever changes.
2074 patches.setdefault((int(record["tract"]), int(record["patch"])), record)
2075 return dict(sorted(patches.items()))
2078def _query_coadd_refs(butler, coadd_dataset, band, keys, skymap_name=None):
2079 """Resolve one dataset ref per ``(tract, patch)`` key, where one exists.
2081 A single constrained query both resolves the ``skymap`` dimension --
2082 which a (visit, detector) data id doesn't carry -- and reports which
2083 patches are absent from the butler's collections.
2085 Parameters
2086 ----------
2087 butler : `lsst.daf.butler.Butler`
2088 Butler to query; its default collections are used.
2089 coadd_dataset : `str`
2090 Dataset name of the coadds (e.g. ``"template_coadd"``).
2091 band : `str`
2092 Band to restrict the query to.
2093 keys : `collections.abc.Container` [`tuple` [`int`, `int`]]
2094 The ``(tract, patch)`` pairs to look for.
2095 skymap_name : `str`, optional
2096 Skymap to restrict the query to. Only needed when the coadds
2097 match more than one skymap.
2099 Returns
2100 -------
2101 refs : `dict` [`tuple` [`int`, `int`], `lsst.daf.butler.DatasetRef`]
2102 Refs of the coadds that exist, keyed on ``(tract, patch)``.
2103 """
2104 where = "band = :band AND tract IN (:tracts) AND patch IN (:patches)"
2105 bind = {"band": band,
2106 "tracts": sorted({tract for tract, _ in keys}),
2107 "patches": sorted({patch for _, patch in keys})}
2108 if skymap_name is not None:
2109 where += " AND skymap = :skymap"
2110 bind["skymap"] = skymap_name
2112 refs = {}
2113 skymaps = set()
2114 for ref in butler.query_datasets(coadd_dataset, where=where, bind=bind,
2115 limit=None, explain=False):
2116 key = (ref.dataId["tract"], ref.dataId["patch"])
2117 # tract and patch are constrained independently, so the query
2118 # returns their cross product; drop the pairs we didn't ask for.
2119 if key not in keys:
2120 continue
2121 skymaps.add(ref.dataId["skymap"])
2122 refs[key] = ref
2123 if len(skymaps) > 1:
2124 raise ValueError(f"{coadd_dataset} matched more than one skymap "
2125 f"({sorted(skymaps)}); pass skymap_name to pick one.")
2126 return refs
2129def _coadd_geometry(butler, ref):
2130 """Read the bounding box and WCS of a coadd without its pixels.
2132 Parameters
2133 ----------
2134 butler : `lsst.daf.butler.Butler`
2135 Butler to read from.
2136 ref : `lsst.daf.butler.DatasetRef`
2137 Reference to the coadd, stored as an afw exposure or an
2138 `lsst.images.cells.CellCoadd`.
2140 Returns
2141 -------
2142 bbox : `lsst.geom.Box2I`
2143 Bounding box of the coadd.
2144 wcs : `lsst.afw.geom.SkyWcs`
2145 WCS of the coadd.
2146 """
2147 components = ref.datasetType.storageClass.allComponents()
2148 wcs_component = "wcs" if "wcs" in components else "sky_projection"
2149 wcs = butler.get(ref.makeComponentRef(wcs_component), storageClass="Wcs")
2150 bbox = butler.get(ref.makeComponentRef("bbox"), storageClass="Box2I")
2151 return bbox, wcs
2154def _template_patch_outlines(butler, template_dataset, data_id, wcs, *,
2155 coadd_dataset, band, skymap_name=None):
2156 """Outline and coadd of each patch recorded as a template input.
2158 A template stored as an `lsst.images.DifferenceImage` records the
2159 butler ID of each coadd it used, in its ``templates`` component. A
2160 legacy template records only the tract and patch, in its
2161 ``coaddInputs`` component, so its coadds are looked up by data ID.
2163 Parameters
2164 ----------
2165 butler : `lsst.daf.butler.Butler`
2166 Butler to load from.
2167 template_dataset : `str`
2168 Dataset name of the template (e.g. ``"template_detector"``).
2169 data_id : `dict`
2170 Data id of the template.
2171 wcs : `lsst.images.SkyProjection`
2172 Sky projection of the difference image the outlines are drawn on.
2173 coadd_dataset : `str`
2174 Dataset name of the coadds; used only for a legacy template.
2175 band : `str`
2176 Band of the coadds; used only for a legacy template.
2177 skymap_name : `str`, optional
2178 Skymap of the coadds; used only for a legacy template.
2180 Returns
2181 -------
2182 outlines : `dict` [`tuple` [`int`, `int`], `list` [`tuple`]]
2183 ``(x, y)`` corners of each patch in difference-image pixels, keyed
2184 on ``(tract, patch)`` and sorted. For a coadd missing from the
2185 butler, the outline of an `lsst.images` template covers only the
2186 part of the patch that overlaps the difference image.
2187 refs : `dict` [`tuple` [`int`, `int`], `lsst.daf.butler.DatasetRef`]
2188 Refs of the coadds that exist, keyed on ``(tract, patch)``.
2189 """
2190 components = butler.get_dataset_type(template_dataset).storageClass.allComponents()
2191 if "templates" not in components: 2191 ↛ 2192line 2191 didn't jump to line 2192 because the condition on line 2191 was never true
2192 patches = _coadd_input_patches(butler, template_dataset, data_id)
2193 outlines = {key: _project_bbox_corners(record.getBBox(), record.getWcs(), wcs)
2194 for key, record in patches.items()}
2195 refs = (_query_coadd_refs(butler, coadd_dataset, band, patches, skymap_name)
2196 if patches else {})
2197 return outlines, refs
2199 outlines = {}
2200 refs = {}
2201 templates = butler.get(f"{template_dataset}.templates", data_id)
2202 for info in sorted(templates, key=lambda info: (info.tract, info.patch)):
2203 key = (info.tract, info.patch)
2204 ref = butler.get_dataset(info.dataset_id)
2205 if ref is None:
2206 # The bounds are already in difference-image pixels; drop the
2207 # closing vertex, which repeats the first.
2208 corners = list(zip(info.bounds.x_vertices.tolist(), info.bounds.y_vertices.tolist()))
2209 outlines[key] = corners[:-1] if corners[0] == corners[-1] else corners
2210 else:
2211 refs[key] = ref
2212 outlines[key] = _project_bbox_corners(*_coadd_geometry(butler, ref), wcs)
2213 return outlines, refs
2216def _project_bbox_corners(bbox, from_wcs, to_wcs):
2217 """Corners of a bounding box mapped from one pixel grid to another.
2219 ``bbox`` (an `lsst.images.Box` or `lsst.geom.Box2I`) is in
2220 ``from_wcs``'s pixel system; the returned ``(x, y)`` corners, taken at
2221 the outer pixel edges, are in ``to_wcs``'s parent pixel system. Each
2222 transform is an `lsst.images.SkyProjection` or `lsst.afw.geom.SkyWcs`.
2223 """
2224 x_min, x_max, y_min, y_max = _box_edges(bbox)
2225 ra, dec = _pixel_to_sky(from_wcs, [x_min, x_max, x_max, x_min], [y_min, y_min, y_max, y_max])
2226 xs, ys = _sky_to_pixel(to_wcs, ra, dec)
2227 return [(float(x), float(y)) for x, y in zip(xs, ys)]
2230def _rect_from_bbox(bbox):
2231 """``(xmin, xmax, ymin, ymax)`` clip rectangle of a bounding box.
2233 Taken at the box's outer pixel edges so that clipped areas are
2234 comparable with the box's area; the integer limits would be a half
2235 pixel short on each side.
2236 """
2237 return _box_edges(bbox)
2240def _draw_outline_on_current_frame(afw_display, corners, ctype, *, label=None,
2241 label_size=1.5, clip_rect=None):
2242 """Draw a closed polyline through ``corners`` on the active frame.
2244 ``clip_rect`` is the displayed image's ``(xmin, xmax, ymin, ymax)``;
2245 when supplied, the label is anchored at the centroid of the outline's
2246 *visible* portion so it stays on-screen for outlines that mostly fall
2247 outside the frame.
2248 """
2249 afw_display.line(list(corners) + [corners[0]], ctype=ctype)
2250 if label is None:
2251 return
2252 visible = _clip_polygon_to_rect(corners, *clip_rect) if clip_rect else corners
2253 if not visible:
2254 visible = corners
2255 x = sum(p[0] for p in visible)/len(visible)
2256 y = sum(p[1] for p in visible)/len(visible)
2257 afw_display.dot(label, x, y, size=label_size, ctype=ctype)
2260def _subset_coadd_to_outline(coadd, corners, margin):
2261 """Trim a coadd to ``margin`` pixels around a projected outline.
2263 Returns the untrimmed coadd if the outline misses it entirely.
2264 """
2265 box = lsst.geom.Box2D()
2266 for x, y in corners:
2267 box.include(lsst.geom.Point2D(x, y))
2268 box.grow(margin)
2269 bbox = lsst.geom.Box2I(box)
2270 bbox.clip(coadd.getBBox())
2271 if bbox.isEmpty():
2272 return coadd
2273 return coadd[bbox]
2276def _print_coadd_coverage_legend(header, coverages, indent=" "):
2277 """Print one row per contributing patch: frame, id, color, coverage."""
2278 print(header)
2279 print(f"{indent}{'frame':>5s} {'tract':>6s} {'patch':>5s} "
2280 f"{'color':<8s} {'overlap':>7s} coadd")
2281 for cov in coverages:
2282 frame = "--" if cov.frame is None else str(cov.frame)
2283 status = "found" if cov.ref is not None else "MISSING from collections"
2284 print(f"{indent}{frame:>5s} {cov.tract:6d} {cov.patch:5d} "
2285 f"{cov.color:<8s} {100*cov.overlap:6.1f}% {status}")
2288def display_coadd_coverage(butler, visit, detector, backend="firefly", *,
2289 patch_extent="full",
2290 patch_margin=100,
2291 show_diffim_outline=True,
2292 show_patch_outlines=True,
2293 diffim_ctype="red",
2294 label_size=1.5,
2295 mask_transparency=80,
2296 strip_metadata=True,
2297 align="Standard",
2298 skymap_name=None,
2299 image_datasets=_IMAGE_DATASETS,
2300 coadd_dataset="template_coadd",
2301 dry_run=False):
2302 """Display a difference image alongside every coadd patch that went
2303 into its template.
2305 Frame 0 shows the difference image, with the outline of each
2306 contributing patch drawn and labeled ``tract,patch``. Frames 1..N show
2307 the coadd patches themselves, one per frame in ``(tract, patch)``
2308 order, each with the difference image's footprint outlined on it.
2310 The contributing patches are read from the template's ``templates``
2311 component (an `lsst.images.DifferenceImage`) or ``coaddInputs``
2312 component (a legacy exposure), which `~lsst.ip.diffim.GetTemplateTask`
2313 fills with one record per patch that supplied valid pixels. That is a
2314 stricter set than a skymap
2315 lookup would give: patches that overlap the detector geometrically but
2316 contributed nothing were already dropped when the template was built.
2318 Parameters
2319 ----------
2320 butler : `lsst.daf.butler.Butler`
2321 Butler to load data from.
2322 visit, detector : `int`
2323 Visit and detector ids to load data for.
2324 backend : `str`, optional
2325 afw display backend (typically "firefly" or "ds9").
2326 patch_extent : {"full", "overlap"}, optional
2327 How much of each coadd to display. ``"full"`` (default) shows the
2328 whole patch, which is what puts the detector's footprint in
2329 context but sends a full-size coadd to the backend per frame;
2330 ``"overlap"`` trims each patch to the difference image's footprint
2331 plus ``patch_margin``, which is much faster to display but makes
2332 every frame look alike.
2333 patch_margin : `int`, optional
2334 Pixels of coadd to keep around the difference image's footprint
2335 when ``patch_extent="overlap"``. Ignored for ``"full"``.
2336 show_diffim_outline : `bool`, optional
2337 If True, outline the difference image's footprint on each coadd
2338 frame, labeled ``visit,detector``.
2339 show_patch_outlines : `bool`, optional
2340 If True, outline each contributing patch on the difference-image
2341 frame, labeled ``tract,patch``. Outlines are drawn for missing
2342 coadds too, from the template's records; for an `lsst.images`
2343 template such an outline covers only the part of the patch that
2344 overlaps the difference image.
2345 diffim_ctype : `str`, optional
2346 Display color for the difference-image outline. The patch outlines
2347 cycle through a fixed palette instead, so that a patch's color on
2348 frame 0 identifies its own frame in the printed legend.
2349 label_size : `float`, optional
2350 Text size for the outline labels.
2351 mask_transparency : `int` or `None`, optional
2352 Mask-plane transparency forwarded to the display (0 = opaque,
2353 100 = fully transparent). Pass ``None`` to leave the backend's
2354 current setting untouched.
2355 strip_metadata : `bool`, optional
2356 Drop ``LTV1``/``LTV2`` keywords from each exposure's metadata
2357 before sending to the backend. Needed for ds9 to align frames.
2358 align : `str` or `None`, optional
2359 ``match_type`` passed to ``alignImages``. Defaults to
2360 ``"Standard"`` (align by WCS) rather than the pixel matching
2361 `display_images` uses, because these frames do not share a pixel
2362 grid. Pass None to leave the frames unaligned.
2363 skymap_name : `str`, optional
2364 Skymap the coadds live in. Only needed when the contributing
2365 tracts and patches match coadds in more than one skymap, which is
2366 an error otherwise.
2367 image_datasets : `dict` [`str`, `str`], optional
2368 Mapping from image-type key (``"science"``, ``"template"``,
2369 ``"difference"``) to butler dataset name. Only the ``"template"``
2370 and ``"difference"`` entries are used here; only the template's
2371 record of its inputs is read, so its pixels are never loaded.
2372 coadd_dataset : `str`, optional
2373 Dataset name of the coadds the template was built from. The
2374 default matches what ``ApPipe.yaml`` binds to the template task's
2375 ``coaddExposures`` input. Ignored for an `lsst.images` template,
2376 which records the butler ID of each coadd.
2377 dry_run : `bool`, optional
2378 If True, work out which patches contributed and print the legend,
2379 but skip constructing the afw display and loading any pixels.
2380 Useful for checking a template's provenance without a viewer.
2381 Default False.
2382 """
2383 if patch_extent not in ("full", "overlap"):
2384 raise ValueError("patch_extent must be 'full' or 'overlap', "
2385 f"got {patch_extent!r}")
2386 data_id = {"visit": visit, "detector": detector}
2387 header = f"visit={visit}, detector={detector} -- "
2389 diffim = _get_image(butler, image_datasets, "difference", data_id)
2390 bbox = diffim.bbox
2391 clip_rect = _rect_from_bbox(bbox)
2392 diffim_area = bbox.area
2393 band = diffim.band
2394 outlines, refs = _template_patch_outlines(
2395 butler, image_datasets["template"], data_id, diffim.sky_projection,
2396 coadd_dataset=coadd_dataset, band=band, skymap_name=skymap_name)
2397 if not outlines:
2398 print(f"{header}template records no coadd inputs; showing the "
2399 "difference image only.")
2401 coverages = []
2402 frame = 1
2403 for i, ((tract, patch), corners) in enumerate(outlines.items()):
2404 overlap = _polygon_area(_clip_polygon_to_rect(corners, *clip_rect))
2405 ref = refs.get((tract, patch))
2406 coverages.append(_PatchCoverage(
2407 tract=tract, patch=patch, ref=ref,
2408 corners_xy=corners, overlap=overlap/diffim_area,
2409 color=_FLAG_PALETTE[i % len(_FLAG_PALETTE)],
2410 frame=None if ref is None else frame))
2411 if ref is not None:
2412 frame += 1
2414 if coverages:
2415 n_tracts = len({cov.tract for cov in coverages})
2416 # Overlaps sum to more than 100%: adjacent patches share a border
2417 # region by construction.
2418 _print_coadd_coverage_legend(
2419 f"{header}{len(coverages)} coadd patches from {n_tracts} tract(s), "
2420 f"band={band}:", coverages)
2422 if dry_run:
2423 return
2425 afw_display = lsst.afw.display.Display(backend=backend)
2426 if mask_transparency is not None:
2427 afw_display.setMaskTransparency(mask_transparency)
2429 afw_display.frame = 0
2430 # Wipe any markers left over from a previous call — `image()`
2431 # only replaces the pixel data, region overlays persist otherwise.
2432 _erase_current_frame_regions(afw_display)
2433 _show_image(afw_display, diffim, "difference", strip_metadata)
2434 if show_patch_outlines:
2435 with afw_display.Buffering():
2436 for cov in coverages:
2437 _draw_outline_on_current_frame(
2438 afw_display, cov.corners_xy, cov.color,
2439 label=f"{cov.tract},{cov.patch}", label_size=label_size,
2440 clip_rect=clip_rect)
2442 for cov in coverages:
2443 if cov.ref is None:
2444 continue
2445 # A coadd is only displayed, so it is read in the legacy form that
2446 # afw displays take.
2447 coadd = _legacy_for_display(butler.get(cov.ref))
2448 # Project after loading rather than reusing the record's wcs, so
2449 # the outline is tied to the pixels actually on screen.
2450 corners = _project_bbox_corners(bbox, diffim.sky_projection, coadd.wcs)
2451 if patch_extent == "overlap":
2452 coadd = _subset_coadd_to_outline(coadd, corners, patch_margin)
2453 afw_display.frame = cov.frame
2454 _erase_current_frame_regions(afw_display)
2455 _show_image(afw_display, coadd, f"{cov.tract},{cov.patch}", strip_metadata)
2456 if show_diffim_outline:
2457 _draw_outline_on_current_frame(
2458 afw_display, corners, diffim_ctype,
2459 label=f"{visit},{detector}", label_size=label_size,
2460 clip_rect=_rect_from_bbox(coadd.getBBox()))
2462 if align is not None:
2463 try:
2464 afw_display.alignImages(match_type=align)
2465 except NotImplementedError:
2466 print(f"WARNING: cannot automatically align and lock images with backend={backend!r}.")
2469def _subset_catalog(catalog, indices):
2470 """Build a `~lsst.afw.table.SourceCatalog` of the rows at `indices`.
2472 The subset shares the input's table and appends the existing records
2473 (shallow), so each row keeps its attached `Footprint` rather than a
2474 copy. Used to split a catalog into one sub-catalog per overlay color.
2475 """
2476 subset = lsst.afw.table.SourceCatalog(catalog.table)
2477 for i in indices:
2478 subset.append(catalog[i])
2479 return subset
2482def _footprint_adjacency(bboxes):
2483 """Adjacency lists for footprints whose bounding boxes overlap.
2485 Two footprints are adjacent when their bounding boxes overlap (share
2486 at least one pixel). Bounding-box overlap can slightly over-count
2487 versus true pixel touching, which only makes the coloring more
2488 conservative (never assigns the same color to two touching footprints).
2490 Parameters
2491 ----------
2492 bboxes : `list` [`lsst.geom.Box2I`]
2493 Footprint bounding boxes, in the order the catalog iterates.
2495 Returns
2496 -------
2497 adjacency : `list` [`set` [`int`]]
2498 ``adjacency[i]`` is the set of indices whose footprints touch
2499 footprint ``i``.
2500 """
2501 n = len(bboxes)
2502 adjacency = [set() for _ in range(n)]
2503 for i in range(n):
2504 for j in range(i + 1, n):
2505 if bboxes[i].overlaps(bboxes[j]):
2506 adjacency[i].add(j)
2507 adjacency[j].add(i)
2508 return adjacency
2511def _greedy_color(adjacency, n_colors, rng=None):
2512 """Greedily color a graph so touching nodes differ, using `n_colors`.
2514 Nodes are colored highest-degree first (Welsh--Powell), and each node
2515 takes a color chosen *at random* from those not already used by an
2516 adjacent node. Randomizing the choice -- rather than always taking the
2517 lowest free index -- spreads the coloring across the whole palette and
2518 varies it between runs, while still guaranteeing adjacent footprints
2519 differ. Footprint adjacency graphs are essentially planar, so with 12
2520 colors available a conflict-free coloring is found in practice. In the
2521 pathological case where a node's neighbors already occupy all
2522 `n_colors`, it falls back to a color least used among those neighbors
2523 (ties broken randomly) rather than failing.
2525 Parameters
2526 ----------
2527 adjacency : `list` [`set` [`int`]]
2528 Adjacency lists from `_footprint_adjacency`.
2529 n_colors : `int`
2530 Number of available colors (palette length).
2531 rng : `random.Random`, optional
2532 Random source, for reproducible colorings in tests. Defaults to a
2533 fresh unseeded `random.Random`.
2535 Returns
2536 -------
2537 colors : `list` [`int`]
2538 ``colors[i]`` is the palette index assigned to node ``i``.
2539 """
2540 if rng is None:
2541 rng = random.Random()
2542 n = len(adjacency)
2543 colors = [-1] * n
2544 order = sorted(range(n), key=lambda i: len(adjacency[i]), reverse=True)
2545 for node in order:
2546 used = {colors[nbr] for nbr in adjacency[node] if colors[nbr] >= 0}
2547 available = [c for c in range(n_colors) if c not in used]
2548 if available:
2549 chosen = rng.choice(available)
2550 else:
2551 # Every color is taken by a neighbor (needs >n_colors mutually
2552 # touching footprints -- essentially never). Reuse a color that
2553 # appears least among the neighbors, breaking ties randomly.
2554 counts = [0] * n_colors
2555 for nbr in adjacency[node]:
2556 if colors[nbr] >= 0:
2557 counts[colors[nbr]] += 1
2558 fewest = min(counts)
2559 chosen = rng.choice(
2560 [c for c in range(n_colors) if counts[c] == fewest])
2561 colors[node] = chosen
2562 return colors
2565def _overlay_footprint_layers(afw_display, buckets, *, style, layer_prefix):
2566 """Overlay footprint buckets as Firefly layers on the current frame.
2568 Firefly's ``overlayFootprints`` takes a single color per call, so a
2569 multi-color overlay is drawn as one layer per bucket. ``layer_prefix``
2570 is prepended to each layer/title string so overlays drawn on different
2571 frames (e.g. the three frames of `display_images`) get distinct layer
2572 ids; the backend appends the frame number itself. The per-bucket suffix
2573 is the bucket's position in ``buckets``, so re-running with the same
2574 number of buckets overwrites the layers in place. ``overlayFootprints``
2575 is a Firefly-impl method reached via `Display`'s attribute delegation,
2576 the same way `display_images` calls ``alignImages``.
2578 Parameters
2579 ----------
2580 afw_display : `lsst.afw.display.Display`
2581 Display with the target frame already selected.
2582 buckets : `list` [`tuple`]
2583 ``(catalog, color, label)`` triples; each ``catalog`` is a
2584 footprint-bearing `~lsst.afw.table.SourceCatalog` drawn in
2585 ``color``. Empty buckets are skipped.
2586 style : {"outline", "fill"}
2587 Footprint rendering style.
2588 layer_prefix : `str`
2589 Prefix for the Firefly layer/title strings; must be unique per
2590 frame to avoid layers on different frames colliding.
2592 Returns
2593 -------
2594 summary : `list` [`tuple`]
2595 ``(label, color, count)`` per drawn bucket, for the legend.
2596 """
2597 summary = []
2598 for i, (catalog, color, label) in enumerate(buckets):
2599 if catalog is None or len(catalog) == 0:
2600 continue
2601 layer = f"{layer_prefix} c{i} "
2602 afw_display.overlayFootprints(catalog, color=color, style=style,
2603 layerString=layer, titleString=layer)
2604 summary.append((label, color, len(catalog)))
2605 return summary
2608def display_footprints(butler=None, visit=None, detector=None,
2609 backend="firefly", *,
2610 exposure=None, catalog=None,
2611 image_type="difference",
2612 catalog_dataset="dia_source_unfiltered",
2613 style="outline",
2614 palette=_OBJECT_PALETTE,
2615 frame=0,
2616 mask_transparency=80,
2617 strip_metadata=True,
2618 image_datasets=_IMAGE_DATASETS):
2619 """Overlay diaSource footprints on an exposure in Firefly, color-cycled
2620 so that touching footprints get distinct colors.
2622 The footprints come from an afw `~lsst.afw.table.SourceCatalog` (the
2623 diffim detection output, which still carries per-source `Footprint`\\ s
2624 -- the transformed and APDB diaSource tables are DataFrames with the
2625 footprints stripped). Supply the data one of two ways:
2627 * pass ``butler`` plus ``visit`` and ``detector`` to load the
2628 exposure (``image_datasets[image_type]``) and catalog
2629 (``catalog_dataset``) from the butler; or
2630 * pass ``exposure`` and ``catalog`` directly, skipping the butler.
2632 The footprints are then drawn on a single Firefly frame using the
2633 backend's native footprint overlay.
2635 Each footprint is assigned one of the `palette` colors by greedy
2636 graph coloring over a bounding-box-touch adjacency graph, so no two
2637 touching footprints share a color (see `_footprint_adjacency` and
2638 `_greedy_color`). The color chosen for each footprint is randomized
2639 among those its neighbors are not using, so the palette is spread
2640 across the frame and re-running produces a different coloring. Because
2641 Firefly's ``overlayFootprints`` takes a single color per call, the
2642 catalog is split into one sub-catalog per color and each is overlaid as
2643 its own Firefly layer.
2645 Re-running on the same frame overwrites each color layer in place.
2646 Color layers left over from a previous run that used *more* colors are
2647 not cleared automatically.
2649 Parameters
2650 ----------
2651 butler : `lsst.daf.butler.Butler`, optional
2652 Butler to load the exposure and catalog from. Required (with
2653 ``visit`` and ``detector``) unless ``exposure`` and ``catalog`` are
2654 given directly.
2655 visit, detector : `int`, optional
2656 Visit and detector ids to load data for. Required with ``butler``.
2657 backend : `str`, optional
2658 afw display backend. Only ``"firefly"`` is supported, since the
2659 overlay uses Firefly's native footprint rendering.
2660 exposure : `lsst.images.VisitImage` or `lsst.afw.image.Exposure`, optional
2661 Image to draw on, supplied directly instead of via the butler.
2662 Must be given together with ``catalog``; when set, ``butler``,
2663 ``visit``, ``detector``, ``catalog_dataset``, ``image_type``, and
2664 ``image_datasets`` are all ignored.
2665 catalog : `lsst.afw.table.SourceCatalog`, optional
2666 Footprint-bearing source catalog, supplied directly instead of via
2667 the butler. Must be given together with ``exposure``.
2668 image_type : {"science", "template", "difference"}, optional
2669 Which image to display the footprints on (butler mode only).
2670 Default ``"difference"``.
2671 catalog_dataset : `str`, optional
2672 Butler dataset of the footprint-bearing afw source catalog (butler
2673 mode only). Default ``"dia_source_unfiltered"`` (the pre-filter
2674 detection catalog, which still carries footprints; the
2675 transformed/standardized diaSource tables have them stripped).
2676 style : {"outline", "fill"}, optional
2677 Footprint rendering style. ``"outline"`` (default) keeps the
2678 color coding legible where footprints overlap; ``"fill"`` shades
2679 the interior.
2680 palette : sequence of `str`, optional
2681 Colors cycled across footprints. Defaults to the 12-color
2682 ``_OBJECT_PALETTE`` also used by the cutout plotters.
2683 frame : `int`, optional
2684 Display frame to draw the image and footprints in. Default ``0``.
2685 mask_transparency : `int` or `None`, optional
2686 Mask-plane transparency forwarded to the display (0 = opaque,
2687 100 = fully transparent). Pass ``None`` to leave it untouched.
2688 strip_metadata : `bool`, optional
2689 Drop ``LTV1``/``LTV2`` keywords from the exposure metadata before
2690 sending to the backend.
2691 image_datasets : `dict` [`str`, `str`], optional
2692 Mapping from image-type key to butler dataset name.
2693 """
2694 if backend != "firefly":
2695 raise ValueError(
2696 f"display_footprints only supports the 'firefly' backend "
2697 f"(needs Firefly's native footprint overlay); got {backend!r}")
2698 if style not in ("outline", "fill"):
2699 raise ValueError(f"style must be 'outline' or 'fill', got {style!r}")
2701 # Two input modes: direct (exposure + catalog) or butler-loaded.
2702 direct = exposure is not None or catalog is not None
2703 if direct:
2704 if exposure is None or catalog is None:
2705 raise ValueError(
2706 "supply BOTH exposure and catalog to draw directly")
2707 title = "footprints"
2708 location = ""
2709 else:
2710 if butler is None or visit is None or detector is None:
2711 raise ValueError(
2712 "supply either (butler, visit, detector) or "
2713 "(exposure, catalog)")
2714 if image_type not in image_datasets:
2715 raise ValueError(
2716 f"image_type must be one of {sorted(image_datasets)}, "
2717 f"got {image_type!r}")
2718 data_id = {"visit": visit, "detector": detector}
2719 exposure = _get_image(butler, image_datasets, image_type, data_id)
2720 catalog = butler.get(catalog_dataset, data_id)
2721 title = f"{image_type} footprints"
2722 location = f"visit={visit}, detector={detector}: "
2724 if not isinstance(catalog, lsst.afw.table.SourceCatalog):
2725 raise TypeError(
2726 f"catalog is a {type(catalog).__name__}, not an afw "
2727 "SourceCatalog. Footprints are only carried by the afw detection "
2728 "catalog (storageClass 'SourceCatalog'); the "
2729 "transformed/standardized diaSource tables (DataFrame or "
2730 "ArrowAstropy, e.g. 'dia_source_detector') have them stripped.")
2731 if len(catalog) > 0 and catalog[0].getFootprint() is None:
2732 raise ValueError(
2733 "catalog is an afw SourceCatalog but its records have no "
2734 "Footprint attached, so there is nothing to draw.")
2736 # Color the footprints so touching ones differ, then group indices by
2737 # assigned color for the per-color Firefly overlay calls below.
2738 bboxes = [record.getFootprint().getBBox() for record in catalog]
2739 color_indices = _greedy_color(_footprint_adjacency(bboxes),
2740 len(palette))
2741 groups = {}
2742 for i, c in enumerate(color_indices):
2743 groups.setdefault(c, []).append(i)
2745 afw_display = lsst.afw.display.Display(backend=backend)
2746 if mask_transparency is not None:
2747 afw_display.setMaskTransparency(mask_transparency)
2748 afw_display.frame = frame
2749 # image() only replaces pixel data; wipe stale region markers first.
2750 _erase_current_frame_regions(afw_display)
2751 _show_image(afw_display, exposure, title, strip_metadata)
2753 print(f"{location}{len(catalog)} footprints in {len(groups)} colors")
2754 buckets = [(_subset_catalog(catalog, groups[c]), palette[c], f"c{c}")
2755 for c in sorted(groups)]
2756 _overlay_footprint_layers(afw_display, buckets, style=style,
2757 layer_prefix="footprints")
2760def extract_timestamped_messages(log: str | dict[str, Any]) -> str:
2761 """Extract ``records[*].(asctime, message)`` from an LSST-style JSON
2762 log and format them one per line, as::
2764 2026-02-25T04:15:35.092108Z Preparing execution...
2766 Parameters
2767 ----------
2768 log : `str` or `dict` [`str`, `Any`]
2769 Either the JSON text or an already-parsed dict.
2771 Returns
2772 -------
2773 messages : `str`
2774 The joined ``asctime message`` lines.
2775 """
2776 if isinstance(log, str):
2777 s = log.strip()
2779 # Handle the case where the *JSON itself* is wrapped in quotes, like:
2780 # '"{...}"' or "'{...}'"
2781 if (len(s) >= 2) and (s[0] == s[-1]) and s[0] in ("'", '"'):
2782 s = s[1:-1]
2784 try:
2785 obj = json.loads(s)
2786 except json.JSONDecodeError:
2787 # One more attempt: sometimes a quoted-JSON string is itself
2788 # JSON-encoded e.g. "\"{...}\""
2789 obj = json.loads(json.loads(s))
2790 else:
2791 obj = log
2793 records = obj.get("records", [])
2794 if not isinstance(records, list):
2795 raise TypeError("Expected obj['records'] to be a list.")
2797 rows: list[tuple[datetime, str, str]] = []
2798 for rec in records:
2799 if not isinstance(rec, dict):
2800 continue
2801 ts = rec.get("asctime")
2802 msg = rec.get("message")
2803 if not ts or msg is None:
2804 continue
2806 # Parse ISO-8601 with trailing "Z"
2807 dt = datetime.fromisoformat(ts.replace("Z", "+00:00")).astimezone(timezone.utc)
2808 rows.append((dt, ts, str(msg)))
2810 return "\n".join(f"{ts} {msg}" for _, ts, msg in rows)