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