Coverage for python/lsst/analysis/ap/skymapOverlay.py: 4%
175 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-22 10:41 +0000
1# This file is part of analysis_ap.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22"""Backend-agnostic skymap tract/patch overlays.
24The geometry of projecting overlapping tract/patch boundaries into a
25display's pixel coordinate system is shared by two very different
26renderers:
28- `draw_skymap_outlines_mpl`, which draws onto a Matplotlib axis whose
29 pixel<->sky relationship is only known approximately (e.g. the
30 spatially-sampled-metrics panels, which carry sampled ``(x, y)`` and
31 ``(coord_ra, coord_dec)`` but no WCS). Use `make_affine_sky_to_xy` to
32 build the required ``sky_to_xy`` callable from those samples.
33- `draw_skymap_outlines_afw`, which draws onto an `lsst.afw.display`
34 frame showing an exposure with a real WCS, so the sky<->pixel mapping
35 is exact.
37Both renderers share `compute_tract_patch_outlines`, which does the
38findTractPatchList lookup, projects every patch corner through a caller
39supplied ``sky_to_xy`` map, and ranks the overlapping tracts by how much
40of the display they cover.
41"""
43from __future__ import annotations
45__all__ = ["make_affine_sky_to_xy", "compute_tract_patch_outlines",
46 "draw_skymap_outlines_mpl", "draw_skymap_outlines_afw"]
48import numpy as np
51def make_affine_sky_to_xy(ra, dec, x, y):
52 """Build a least-squares affine map from sky to detector pixels.
54 Useful when no WCS is available but matched ``(ra, dec)`` and
55 ``(x, y)`` samples are. For a single detector the affine fit is
56 typically accurate to well under a pixel -- enough for visualization
57 but not for science.
59 Parameters
60 ----------
61 ra, dec : array-like
62 Sky coordinates of the samples, in **radians**.
63 x, y : array-like
64 Detector pixel coordinates of the same samples.
66 Returns
67 -------
68 sky_to_xy : callable
69 Maps an `lsst.geom.SpherePoint` to an ``(x, y)`` tuple of
70 `float` in the same pixel system as the input ``x, y``.
71 """
72 ra = np.asarray(ra)
73 dec = np.asarray(dec)
74 A = np.column_stack([np.ones(ra.size), ra, dec])
75 coef_x, *_ = np.linalg.lstsq(A, np.asarray(x), rcond=None)
76 coef_y, *_ = np.linalg.lstsq(A, np.asarray(y), rcond=None)
78 def sky_to_xy(sphere_point):
79 r = sphere_point.getRa().asRadians()
80 d = sphere_point.getDec().asRadians()
81 return (float(coef_x[0] + coef_x[1]*r + coef_x[2]*d),
82 float(coef_y[0] + coef_y[1]*r + coef_y[2]*d))
84 return sky_to_xy
87def compute_tract_patch_outlines(skymap, sky_to_xy, sky_corners, clip_rect):
88 """Project overlapping tract/patch boundaries into display coordinates.
90 Parameters
91 ----------
92 skymap : `lsst.skymap.BaseSkyMap`
93 Skymap to query for overlapping tracts and patches.
94 sky_to_xy : callable
95 Maps an `lsst.geom.SpherePoint` to an ``(x, y)`` tuple in the
96 display's pixel coordinate system.
97 sky_corners : `list` [`lsst.geom.SpherePoint`]
98 Sky positions spanning the region of interest; passed to
99 ``skymap.findTractPatchList`` to enumerate the tracts/patches that
100 touch it.
101 clip_rect : `tuple` [`float`]
102 ``(xmin, xmax, ymin, ymax)`` display-pixel rectangle. Used only to
103 rank tracts by their visible (clipped) patch area, so that the
104 most-covering tract sorts first; it does not clip the returned
105 geometry.
107 Returns
108 -------
109 outlines : `list` [`dict`]
110 One entry per overlapping tract, sorted by descending visible
111 area (ties broken by tract id for determinism). Each is::
113 {"tract_id": int,
114 "rank": int, # 0-based position in the sorted list
115 "patches": [{"patch_index": int, # sequential index
116 "corners_xy": [(x, y), ...], # 4 corners
117 "center_xy": (x, y)}, # inner-bbox center
118 ...]}
119 """
120 import lsst.geom as geom
122 xmin, xmax, ymin, ymax = clip_rect
123 tract_patch_list = skymap.findTractPatchList(sky_corners)
125 tract_entries = [] # (visible_area, tract_id, patches)
126 for tract_info, patches in tract_patch_list:
127 tract_wcs = tract_info.wcs
128 per_tract_patches = []
129 per_tract_area = 0.0
130 for patch in patches:
131 bbox = patch.getInnerBBox()
132 corners_px = [geom.Point2D(bbox.minX, bbox.minY),
133 geom.Point2D(bbox.maxX, bbox.minY),
134 geom.Point2D(bbox.maxX, bbox.maxY),
135 geom.Point2D(bbox.minX, bbox.maxY)]
136 xy_corners = [sky_to_xy(tract_wcs.pixelToSky(p)) for p in corners_px]
137 center_px = geom.Point2D(0.5*(bbox.minX + bbox.maxX),
138 0.5*(bbox.minY + bbox.maxY))
139 center_xy = sky_to_xy(tract_wcs.pixelToSky(center_px))
141 clipped = _clip_polygon_to_rect(xy_corners, xmin, xmax, ymin, ymax)
142 per_tract_area += _polygon_area(clipped)
143 per_tract_patches.append({"patch_index": patch.getSequentialIndex(),
144 "corners_xy": xy_corners,
145 "center_xy": center_xy})
146 tract_entries.append((per_tract_area, tract_info.getId(), per_tract_patches))
148 # Sort by descending visible area; break ties by tract id so the order
149 # (and hence the mpl linestyle assignment) is stable across calls.
150 tract_entries.sort(key=lambda e: (-e[0], e[1]))
151 return [{"tract_id": tract_id, "rank": rank, "patches": patches}
152 for rank, (_area, tract_id, patches) in enumerate(tract_entries)]
155def draw_skymap_outlines_mpl(ax, skymap, sky_to_xy, sky_corners, *,
156 label_fontsize=7):
157 """Overlay patch boundaries (with ``tract,patch`` labels) on a panel.
159 Each overlapping patch gets a thin white outline (black-stroked so it
160 reads on any colormap) plus a ``tract,patch`` label placed along the
161 midpoint of its longest visible edge inside the current view.
163 When more than one tract overlaps the panel (e.g. detectors near a
164 tract boundary), patches are distinguished by linestyle: tracts are
165 ranked by their total visible patch area and assigned solid, dashed,
166 dotted, dash-dotted in turn.
168 Parameters
169 ----------
170 ax : `matplotlib.axes.Axes`
171 Axis to draw on. Its current ``xlim``/``ylim`` define the visible
172 region and are restored on return.
173 skymap : `lsst.skymap.BaseSkyMap`
174 Skymap to query for overlapping tracts and patches.
175 sky_to_xy : callable
176 Maps an `lsst.geom.SpherePoint` to an ``(x, y)`` tuple in the
177 axis' data coordinate system (see `make_affine_sky_to_xy`).
178 sky_corners : `list` [`lsst.geom.SpherePoint`]
179 Sky positions spanning the region of interest, used to enumerate
180 the overlapping tracts/patches.
181 label_fontsize : `int` or `float`, optional
182 Font size for the per-patch ``tract,patch`` labels.
183 """
184 import matplotlib.patheffects as pe
186 xlim = ax.get_xlim()
187 ylim = ax.get_ylim()
188 xmin, xmax = xlim
189 ymin, ymax = ylim
190 clip_rect = (xmin, xmax, ymin, ymax)
192 outlines = compute_tract_patch_outlines(skymap, sky_to_xy, sky_corners, clip_rect)
194 # White lines with a thin black stroke read on any colormap background.
195 line_outline = [pe.withStroke(linewidth=2.0, foreground="black")]
196 text_outline = [pe.withStroke(linewidth=1.4, foreground="black")]
197 linestyles = ("-", "--", ":", "-.")
199 for tract in outlines:
200 rank = tract["rank"]
201 tract_id = tract["tract_id"]
202 linestyle = linestyles[rank] if rank < len(linestyles) else "-"
203 patch_kwargs = dict(color="white", linewidth=0.5, alpha=0.85,
204 linestyle=linestyle,
205 path_effects=line_outline, zorder=6)
206 for patch in tract["patches"]:
207 xy_corners = patch["corners_xy"]
208 xs_c = [p[0] for p in xy_corners] + [xy_corners[0][0]]
209 ys_c = [p[1] for p in xy_corners] + [xy_corners[0][1]]
210 ax.plot(xs_c, ys_c, **patch_kwargs)
212 # Place the label at the midpoint of the patch edge with the
213 # longest visible portion within the panel, offset slightly
214 # toward the patch center so the text sits inside.
215 cx, cy = patch["center_xy"]
216 best = None
217 for i in range(4):
218 (x0, y0), (x1, y1) = xy_corners[i], xy_corners[(i + 1) % 4]
219 clipped = _clip_segment_to_rect(x0, y0, x1, y1,
220 xmin, xmax, ymin, ymax)
221 if clipped is None:
222 continue
223 cx0, cy0, cx1, cy1 = clipped
224 length = float(np.hypot(cx1 - cx0, cy1 - cy0))
225 if best is None or length > best[0]:
226 best = (length, cx0, cy0, cx1, cy1)
227 if best is None:
228 continue # entire patch is outside the panel
229 length, cx0, cy0, cx1, cy1 = best
230 mx = 0.5*(cx0 + cx1)
231 my = 0.5*(cy0 + cy1)
232 # Inward offset toward the projected patch center. Use the
233 # smaller of "fixed fraction of edge length" and "fraction of
234 # the midpoint-to-center distance" so the offset never lands
235 # outside the patch on slivers.
236 dx, dy = cx - mx, cy - my
237 d_center = float(np.hypot(dx, dy))
238 if d_center > 0:
239 step = min(0.06*length, 0.4*d_center)
240 mx += dx/d_center * step
241 my += dy/d_center * step
243 # Rotate the text to lie parallel to the visible edge, flipping
244 # to keep it reading right-side-up (angle clamped to [-90, 90]).
245 angle_deg = float(np.degrees(np.arctan2(cy1 - cy0, cx1 - cx0)))
246 if angle_deg > 90.0:
247 angle_deg -= 180.0
248 elif angle_deg < -90.0:
249 angle_deg += 180.0
251 ax.text(mx, my, f"{tract_id},{patch['patch_index']}",
252 ha="center", va="center",
253 rotation=angle_deg, rotation_mode="anchor",
254 color="white", fontsize=label_fontsize,
255 path_effects=text_outline, zorder=7,
256 clip_on=True)
258 ax.set_xlim(xlim)
259 ax.set_ylim(ylim)
262def draw_skymap_outlines_afw(afw_display, skymap, wcs, bbox, *,
263 ctype="green", label_size=1.5, draw_labels=True):
264 """Overlay tract/patch boundaries on the current `afw.display` frame.
266 Each overlapping patch is drawn as a closed polyline in the image's
267 parent pixel coordinates and, optionally, labeled ``tract,patch``
268 near the center of its visible portion.
270 Parameters
271 ----------
272 afw_display : `lsst.afw.display.Display`
273 Display whose current frame already shows the exposure. The caller
274 is responsible for selecting the frame.
275 skymap : `lsst.skymap.BaseSkyMap`
276 Skymap to query for overlapping tracts and patches.
277 wcs : `lsst.afw.geom.SkyWcs`
278 WCS of the displayed exposure (``exposure.wcs``).
279 bbox : `lsst.geom.Box2I` or `lsst.geom.Box2D`
280 Parent bounding box of the displayed exposure.
281 Defines the footprint searched for overlapping patches
282 and the region used to place labels.
283 ctype : `str`, optional
284 Display color for both the outlines and labels.
285 label_size : `float`, optional
286 Text size for the ``tract,patch`` labels.
287 draw_labels : `bool`, optional
288 If False, draw only the outlines.
289 """
290 import lsst.geom as geom
292 def sky_to_xy(sphere_point):
293 p = wcs.skyToPixel(sphere_point)
294 return (p.getX(), p.getY())
296 xmin = bbox.getMinX()
297 xmax = bbox.getMaxX()
298 ymin = bbox.getMinY()
299 ymax = bbox.getMaxY()
300 clip_rect = (xmin, xmax, ymin, ymax)
302 # Sky positions at the four image corners span the footprint for the
303 # tract/patch lookup.
304 corner_px = [geom.Point2D(xmin, ymin), geom.Point2D(xmax, ymin),
305 geom.Point2D(xmax, ymax), geom.Point2D(xmin, ymax)]
306 sky_corners = [wcs.pixelToSky(p) for p in corner_px]
308 outlines = compute_tract_patch_outlines(skymap, sky_to_xy, sky_corners, clip_rect)
310 with afw_display.Buffering():
311 for tract in outlines:
312 tract_id = tract["tract_id"]
313 for patch in tract["patches"]:
314 xy_corners = patch["corners_xy"]
315 # Closed polyline: repeat the first corner.
316 afw_display.line(list(xy_corners) + [xy_corners[0]], ctype=ctype)
317 if not draw_labels:
318 continue
319 # Anchor the label at the centroid of the patch's visible
320 # portion so it doesn't land off-image for patches that
321 # mostly fall outside the detector.
322 clipped = _clip_polygon_to_rect(xy_corners, xmin, xmax, ymin, ymax)
323 if clipped:
324 lx = sum(p[0] for p in clipped)/len(clipped)
325 ly = sum(p[1] for p in clipped)/len(clipped)
326 else:
327 lx, ly = patch["center_xy"]
328 afw_display.dot(f"{tract_id},{patch['patch_index']}",
329 lx, ly, size=label_size, ctype=ctype)
332def _clip_polygon_to_rect(polygon, xmin, xmax, ymin, ymax):
333 """Clip a convex polygon against an axis-aligned rectangle.
335 Parameters
336 ----------
337 polygon : sequence of ``(x, y)`` tuples
338 Vertices of the (convex) input polygon, in order.
339 xmin, xmax, ymin, ymax : `float`
340 The clipping rectangle.
342 Returns
343 -------
344 clipped : `list` of ``(x, y)`` tuples
345 The clipped polygon, or an empty list if the polygon lies
346 entirely outside the rectangle.
347 """
348 # Each clip edge is parameterized by ("axis", value, keep_side)
349 # where keep_side is +1 if "inside" means coordinate >= value,
350 # -1 if "inside" means coordinate <= value.
351 edges = (("x", xmin, +1), ("x", xmax, -1),
352 ("y", ymin, +1), ("y", ymax, -1))
354 def _inside(point, axis, val, sign):
355 coord = point[0] if axis == "x" else point[1]
356 return (coord - val)*sign >= 0.0
358 def _intersect(p1, p2, axis, val):
359 x1, y1 = p1
360 x2, y2 = p2
361 if axis == "x":
362 t = (val - x1)/(x2 - x1)
363 return (val, y1 + t*(y2 - y1))
364 t = (val - y1)/(y2 - y1)
365 return (x1 + t*(x2 - x1), val)
367 output = list(polygon)
368 for axis, val, sign in edges:
369 if not output:
370 return []
371 input_list = output
372 output = []
373 for i in range(len(input_list)):
374 curr = input_list[i]
375 prev = input_list[i - 1]
376 curr_in = _inside(curr, axis, val, sign)
377 prev_in = _inside(prev, axis, val, sign)
378 if curr_in:
379 if not prev_in:
380 output.append(_intersect(prev, curr, axis, val))
381 output.append(curr)
382 elif prev_in:
383 output.append(_intersect(prev, curr, axis, val))
384 return output
387def _polygon_area(polygon):
388 """Area of a polygon via the shoelace formula."""
389 n = len(polygon)
390 if n < 3:
391 return 0.0
392 s = 0.0
393 for i in range(n):
394 x1, y1 = polygon[i]
395 x2, y2 = polygon[(i + 1) % n]
396 s += x1*y2 - x2*y1
397 return abs(s)*0.5
400def _clip_segment_to_rect(x0, y0, x1, y1, xmin, xmax, ymin, ymax):
401 """Liang-Barsky line-segment clipping against an axis-aligned rect.
403 Returns the clipped endpoints ``(x0', y0', x1', y1')`` or ``None`` if
404 the segment lies entirely outside the rectangle.
405 """
406 dx = x1 - x0
407 dy = y1 - y0
408 p = (-dx, dx, -dy, dy)
409 q = (x0 - xmin, xmax - x0, y0 - ymin, ymax - y0)
410 u1, u2 = 0.0, 1.0
411 for pi, qi in zip(p, q):
412 if pi == 0.0:
413 if qi < 0.0:
414 return None
415 else:
416 t = qi/pi
417 if pi < 0.0:
418 if t > u2:
419 return None
420 if t > u1:
421 u1 = t
422 else:
423 if t < u1:
424 return None
425 if t < u2:
426 u2 = t
427 return (x0 + u1*dx, y0 + u1*dy, x0 + u2*dx, y0 + u2*dy)