diff --git a/docs/superpowers/plans/2026-07-03-wcs-axes.md b/docs/superpowers/plans/2026-07-03-wcs-axes.md new file mode 100644 index 0000000..a7c7656 --- /dev/null +++ b/docs/superpowers/plans/2026-07-03-wcs-axes.md @@ -0,0 +1,1071 @@ +# WCS Axes via starlink-pyast Implementation Plan + +> **For agentic workers:** REQUIRED SUB-SKILL: Use superpowers:subagent-driven-development (recommended) or superpowers:executing-plans to implement this plan task-by-task. Steps use checkbox (`- [ ]`) syntax for tracking. + +**Goal:** When `mtv()` displays data with a WCS, all axis furniture (border, ticks, labels, optional grid) is drawn by AST in sky coordinates via starlink-pyast, replacing matplotlib's pixel axes. + +**Architecture:** A new `wcsAxes.py` module converts an afw `SkyWcs` to a pyast `FrameSet` (exact AST-native serialization round-trip, no FITS approximation) and provides a `WcsAxesManager` that draws/redraws AST axes on a matplotlib `Axes` whose data coordinates are LSST pixel coordinates. `DisplayImpl` in `matplotlib.py` creates one manager per axes at `mtv` time and keeps it alive through zoom/pan/erase. + +**Tech Stack:** starlink-pyast (`starlink.Ast`, `starlink.Grf.grf_matplotlib`), astshim, matplotlib, lsst.afw. + +**Spec:** `docs/superpowers/specs/2026-07-03-wcs-axes-design.md` + +## Global Constraints + +- flake8 clean per `setup.cfg`: `max-line-length = 110`, `max-doc-length = 79`; run `flake8 python tests` before every commit and expect no output. +- Top-level imports only (no function-scoped imports) in new code. +- Never push to the remote; commit locally only. +- End every commit message with `Co-Authored-By: Claude Fable 5 `. +- One manager instance per matplotlib `Axes`; AST artists are matplotlib `Line2D`/`Text` objects added to `ax.lines`/`ax.texts`. +- All test/verification commands run inside the EUPS environment (see below). + +### Environment setup (needed once per shell, before any pytest/python command) + +```bash +export LSSTSW=/Users/timj/work/lsstsw +source "$LSSTSW/bin/envconfig" +setup -t b8267 lsst_distrib +cd /Users/timj/work/lsstsw/build/display_matplotlib +setup -k -r . +``` + +Note: shell state does not persist between separate Bash tool calls, so each test command must be prefixed by this whole block in one call. + +### Verified facts (from prototyping; do not re-derive) + +- `wcs.getFrameDict()` serialized via `astshim.Channel(astshim.StringStream())` and read back with `starlink.Ast.Channel` reproduces `pixelToSky` exactly. The FrameSet base frame is LSST 0-based pixel coordinates, identical to the axes data coordinates set by `imshow(extent=...)`. +- `starlink.Ast.Channel`'s source must be an object with an `astsource()` method returning one line per call and `None` at end; a bare callable does not work. +- `grf_matplotlib.Line/Mark/Text` draw in **data coordinates** (no explicit transform), so the Plot's graphics box equals its base box: both are the axes view limits in pixel coordinates, ordered `(xlo, ylo, xhi, yhi)`. +- `grf_matplotlib` needs a renderer: `fig.canvas.get_renderer()` must exist (true for Agg-based canvases, i.e. `pyplot.figure()` under Agg/QtAgg/TkAgg). Exotic canvases without it raise `AttributeError` — the error fallback path covers this. +- `Axes.add_artist` clips artists to the axes patch; AST's exterior labels sit *outside* the patch, so every tracked `Text` artist needs `set_clip_on(False)` after `plot.grid()`. +- AST draws its own title ("ICRS coordinates; gnomonic projection"); pass `DrawTitle=0` so matplotlib's `set_title` is used instead. +- Setting `Format(1)=d` (no precision) disables AST's automatic digit selection and renders every RA label as e.g. "30". Decimal-degree mode must compute precision explicitly (see `_decimalDigits` in Task 2). +- `ax.lines = []` raises `AttributeError` on matplotlib ≥ 3.7 (installed: 3.10.9), i.e. the current `_erase()` is broken; Task 5 fixes it as a side effect. +- `wcs.getPixelScale().asDegrees()` works with no arguments (returns degrees/pixel). +- `afwDisplay.delAllDisplays()` exists for test teardown. + +## File Structure + +- Create: `python/lsst/display/matplotlib/wcsAxes.py` — SkyWcs→FrameSet conversion and `WcsAxesManager` (all AST/Grf knowledge lives here). +- Modify: `python/lsst/display/matplotlib/matplotlib.py` — `DisplayImpl` gains three kwargs and delegates to `WcsAxesManager`; `_erase`/`_close`/`useSexagesimal` become WCS-axes aware. +- Modify: `tests/test_display_matplotlib.py` — all new tests. + +--- + +### Task 1: SkyWcs → pyast FrameSet conversion + +**Files:** +- Create: `python/lsst/display/matplotlib/wcsAxes.py` +- Test: `tests/test_display_matplotlib.py` + +**Interfaces:** +- Consumes: `lsst.afw.geom.SkyWcs` (from afw), `astshim`, `starlink.Ast`. +- Produces: `astFrameSetFromWcs(wcs) -> starlink.Ast.FrameSet` — base frame is LSST 0-based pixel coordinates (domain `PIXELS`), current frame is `SKY` in radians. Task 2 builds on this module; Task 4 imports this function into `matplotlib.py`. + +- [ ] **Step 1: Write the failing test** + +Add to `tests/test_display_matplotlib.py`. Replace the current import block (lines 24-27) with: + +```python +import unittest + +import matplotlib +matplotlib.use("Agg") + +import lsst.utils.tests # noqa: E402 +import lsst.geom # noqa: E402 +import lsst.afw.display as afwDisplay # noqa: E402 +import lsst.afw.geom as afwGeom # noqa: E402 +``` + +(Later tasks add further imports as they become used: `matplotlib.figure`, `matplotlib.text`, `FigureCanvasAgg` in Task 2; `lsst.afw.image` in Task 4. Adding them earlier would fail flake8 F401.) + +Then add after the imports: + +```python +def makeTestWcs(): + """Return a rotated TAN SkyWcs for testing.""" + return afwGeom.makeSkyWcs( + crpix=lsst.geom.Point2D(100, 75), + crval=lsst.geom.SpherePoint(30.0, -45.0, lsst.geom.degrees), + cdMatrix=afwGeom.makeCdMatrix(scale=3.0*lsst.geom.arcseconds, + orientation=30*lsst.geom.degrees), + ) + + +class AstFrameSetTestCase(lsst.utils.tests.TestCase): + """Tests for the SkyWcs to pyast FrameSet conversion.""" + + def testRoundTrip(self): + """The pyast FrameSet must reproduce pixelToSky exactly.""" + from lsst.display.matplotlib.wcsAxes import astFrameSetFromWcs + + wcs = makeTestWcs() + frameSet = astFrameSetFromWcs(wcs) + for x, y in [(0.0, 0.0), (100.0, 75.0), (199.5, 149.5), (-20.0, 300.0)]: + sky = wcs.pixelToSky(x, y) + result = frameSet.tran([[x], [y]]) + self.assertEqual(result[0][0], sky[0].asRadians()) + self.assertEqual(result[1][0], sky[1].asRadians()) +``` + +- [ ] **Step 2: Run test to verify it fails** + +```bash +export LSSTSW=/Users/timj/work/lsstsw +source "$LSSTSW/bin/envconfig" +setup -t b8267 lsst_distrib +cd /Users/timj/work/lsstsw/build/display_matplotlib +setup -k -r . +pytest tests/test_display_matplotlib.py::AstFrameSetTestCase -v +``` + +Expected: FAIL with `ModuleNotFoundError: No module named 'lsst.display.matplotlib.wcsAxes'`. + +- [ ] **Step 3: Write minimal implementation** + +Create `python/lsst/display/matplotlib/wcsAxes.py`: + +```python +# +# LSST Data Management System +# Copyright 2026 LSST Corporation. +# +# This product includes software developed by the +# LSST Project (http://www.lsst.org/). +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the LSST License Statement and +# the GNU General Public License along with this program. If not, +# see . +# + +"""Draw sky coordinate axes on a matplotlib Axes using AST. + +The axes are drawn by `starlink.Ast.Plot` using the matplotlib grf +plugin, so arbitrary AST WCS transforms are supported with no FITS +approximation. +""" + +__all__ = ["astFrameSetFromWcs", "WcsAxesManager"] + +import astshim +import starlink.Ast + + +class _AstStringSource: + """Present AST native serialization text to `starlink.Ast.Channel`. + + `starlink.Ast.Channel` reads via an object with an ``astsource`` + method returning one line per call, and `None` at end of input. + + Parameters + ---------- + text : `str` + AST native serialization of an AST object. + """ + + def __init__(self, text): + self._lines = text.splitlines() + self._pos = 0 + + def astsource(self): + if self._pos >= len(self._lines): + return None + line = self._lines[self._pos] + self._pos += 1 + return line + + +def astFrameSetFromWcs(wcs): + """Convert an afw SkyWcs to a starlink-pyast FrameSet. + + Parameters + ---------- + wcs : `lsst.afw.geom.SkyWcs` + The WCS to convert. + + Returns + ------- + frameSet : `starlink.Ast.FrameSet` + FrameSet whose base frame is LSST 0-based pixel coordinates + (domain ``PIXELS``) and whose current frame is sky coordinates + in radians (domain ``SKY``). + + Notes + ----- + The conversion serializes the WCS's underlying `astshim.FrameDict` + to AST native text and reads it back with pyast, so the result is + exact; no FITS approximation is involved. + """ + stream = astshim.StringStream() + astshim.Channel(stream).write(wcs.getFrameDict()) + channel = starlink.Ast.Channel(_AstStringSource(stream.getSinkData())) + return channel.read() +``` + +(The module docstring mentions `starlink.Ast.Plot`; the Plot itself and its imports arrive with `WcsAxesManager` in Task 2.) + +- [ ] **Step 4: Run test to verify it passes** + +Same command as Step 2. Expected: PASS. + +- [ ] **Step 5: Run flake8 and commit** + +```bash +flake8 python tests +git add python/lsst/display/matplotlib/wcsAxes.py tests/test_display_matplotlib.py +git commit -m "Add exact SkyWcs to pyast FrameSet conversion + +Co-Authored-By: Claude Fable 5 " +``` + +--- + +### Task 2: WcsAxesManager — draw, track, remove + +**Files:** +- Modify: `python/lsst/display/matplotlib/wcsAxes.py` +- Test: `tests/test_display_matplotlib.py` + +**Interfaces:** +- Consumes: `astFrameSetFromWcs` output (Task 1); a matplotlib `Axes` whose data coordinates are LSST pixel coordinates. +- Produces: `WcsAxesManager(axes, frameSet, pixelScaleDeg, *, grid=False, useSexagesimal=False, extraOptions="")` with attribute `artists` (list of matplotlib artists) and methods `draw()` and `remove()`. Task 3 adds `setSexagesimal(useSexagesimal)` and limit-change callbacks; Task 4 constructs it from `DisplayImpl`. + +- [ ] **Step 1: Write the failing tests** + +In `tests/test_display_matplotlib.py`, extend the import block (keeping `matplotlib.use("Agg")` immediately after `import matplotlib`): + +```python +import matplotlib +matplotlib.use("Agg") +import matplotlib.figure # noqa: E402 +import matplotlib.text # noqa: E402 +from matplotlib.backends.backend_agg import FigureCanvasAgg # noqa: E402 +``` + +Then add: + +```python +def makeWcsAxes(**kwargs): + """Create an Agg figure with pixel-coordinate axes and a manager. + + Returns the (axes, manager) pair. + """ + from lsst.display.matplotlib.wcsAxes import WcsAxesManager, astFrameSetFromWcs + + fig = matplotlib.figure.Figure(figsize=(8, 6)) + FigureCanvasAgg(fig) + ax = fig.add_subplot(111) + ax.set_xlim(-0.5, 199.5) + ax.set_ylim(-0.5, 149.5) + wcs = makeTestWcs() + manager = WcsAxesManager(ax, astFrameSetFromWcs(wcs), + wcs.getPixelScale().asDegrees(), **kwargs) + return ax, manager + + +class WcsAxesManagerTestCase(lsst.utils.tests.TestCase): + """Tests for WcsAxesManager drawing and artist lifecycle.""" + + def testDrawCreatesArtists(self): + ax, manager = makeWcsAxes() + self.assertGreater(len(manager.artists), 0) + # All tracked artists really are in the axes. + axesArtists = set(ax.lines) | set(ax.texts) + for artist in manager.artists: + self.assertIn(artist, axesArtists) + # Sky axis labels come from the AST SkyFrame. + texts = [a.get_text() for a in manager.artists + if isinstance(a, matplotlib.text.Text)] + self.assertIn("Right ascension", texts) + self.assertIn("Declination", texts) + + def testFurnitureHidden(self): + ax, manager = makeWcsAxes() + self.assertFalse(ax.xaxis.get_visible()) + self.assertFalse(ax.yaxis.get_visible()) + for spine in ax.spines.values(): + self.assertFalse(spine.get_visible()) + + def testTextNotClipped(self): + ax, manager = makeWcsAxes() + for artist in manager.artists: + if isinstance(artist, matplotlib.text.Text): + self.assertFalse(artist.get_clip_on()) + + def testDecimalLabelsByDefault(self): + ax, manager = makeWcsAxes() + tickLabels = [a.get_text() for a in manager.artists + if isinstance(a, matplotlib.text.Text) + and a.get_text() not in ("Right ascension", "Declination")] + self.assertGreater(len(tickLabels), 0) + for label in tickLabels: + self.assertNotIn(":", label) + # 3" pixels over a 200 pixel span is ~0.17 deg: expect 3 decimals. + self.assertTrue(any("." in label and len(label.split(".")[1]) == 3 + for label in tickLabels)) + + def testSexagesimalLabels(self): + ax, manager = makeWcsAxes(useSexagesimal=True) + tickLabels = [a.get_text() for a in manager.artists + if isinstance(a, matplotlib.text.Text)] + self.assertTrue(any(":" in label for label in tickLabels)) + + def testRemove(self): + ax, manager = makeWcsAxes() + manager.remove() + self.assertEqual(manager.artists, []) + self.assertEqual(len(ax.lines), 0) + self.assertEqual(len(ax.texts), 0) + self.assertTrue(ax.xaxis.get_visible()) + self.assertTrue(ax.yaxis.get_visible()) + for spine in ax.spines.values(): + self.assertTrue(spine.get_visible()) +``` + +- [ ] **Step 2: Run tests to verify they fail** + +```bash +export LSSTSW=/Users/timj/work/lsstsw +source "$LSSTSW/bin/envconfig" +setup -t b8267 lsst_distrib +cd /Users/timj/work/lsstsw/build/display_matplotlib +setup -k -r . +pytest tests/test_display_matplotlib.py::WcsAxesManagerTestCase -v +``` + +Expected: FAIL with `ImportError: cannot import name 'WcsAxesManager'`. + +- [ ] **Step 3: Implement WcsAxesManager** + +In `python/lsst/display/matplotlib/wcsAxes.py`, replace the import block with: + +```python +import math + +import astshim +import matplotlib.text +import starlink.Ast +from starlink.Grf import grf_matplotlib +``` + +Then append: + +```python +class WcsAxesManager: + """Manage AST-drawn sky coordinate axes on a matplotlib Axes. + + All visible axis furniture (border, ticks, labels, optional grid) + is drawn by AST; matplotlib's own axis furniture is hidden while + this manager is active and restored by `remove`. + + The axes' data coordinates must be LSST 0-based pixel coordinates + matching the base frame of ``frameSet``, as produced by + `astFrameSetFromWcs`. + + Parameters + ---------- + axes : `matplotlib.axes.Axes` + The axes to draw on. + frameSet : `starlink.Ast.FrameSet` + Pixel-to-sky transform, base frame in axes data coordinates. + pixelScaleDeg : `float` + Pixel scale in degrees per pixel, used to choose the label + precision for decimal-degree formatting. + grid : `bool`, optional + Draw the full curvilinear coordinate grid across the image? + useSexagesimal : `bool`, optional + Format labels as e.g. HH:MM:SS.s rather than decimal degrees. + extraOptions : `str`, optional + Comma-separated AST Plot attribute settings appended after the + defaults, so they take precedence (e.g. + ``"Colour(grid)=red, Size(TextLab)=12"``). + """ + + def __init__(self, axes, frameSet, pixelScaleDeg, *, + grid=False, useSexagesimal=False, extraOptions=""): + self._axes = axes + self._frameSet = frameSet + self._pixelScaleDeg = pixelScaleDeg + self._grid = grid + self._useSexagesimal = useSexagesimal + self._extraOptions = extraOptions + self.artists = [] + + self._savedVisibility = ( + axes.xaxis.get_visible(), + axes.yaxis.get_visible(), + {name: spine.get_visible() for name, spine in axes.spines.items()}, + ) + axes.xaxis.set_visible(False) + axes.yaxis.set_visible(False) + for spine in axes.spines.values(): + spine.set_visible(False) + + self.draw() + + def _decimalDigits(self): + """Return the number of decimal places for degree labels. + + Chosen so that labels remain distinct across the current view: + two more digits than the order of magnitude of the view's + angular extent. + """ + xlo, xhi = self._axes.get_xlim() + ylo, yhi = self._axes.get_ylim() + span = max(abs(xhi - xlo), abs(yhi - ylo)) + extentDeg = self._pixelScaleDeg*span + if extentDeg <= 0: + return 4 + return int(min(10, max(1, math.ceil(-math.log10(extentDeg)) + 2))) + + def _plotOptions(self): + """Return the AST Plot options string for the current state.""" + options = ["DrawTitle=0"] + options.append("Grid=1" if self._grid else "Grid=0") + if not self._useSexagesimal: + digits = self._decimalDigits() + options.append(f"Format(1)=d.{digits}") + options.append(f"Format(2)=d.{digits}") + if self._extraOptions: + options.append(self._extraOptions) + return ", ".join(options) + + def draw(self): + """Draw (or redraw) the AST axes for the current view limits.""" + self._removeArtists() + + xlo, xhi = self._axes.get_xlim() + ylo, yhi = self._axes.get_ylim() + box = (xlo, ylo, xhi, yhi) + + before = set(self._axes.lines) | set(self._axes.texts) + plot = starlink.Ast.Plot(self._frameSet, box, box, + grf_matplotlib(self._axes), + self._plotOptions()) + plot.grid() + self.artists = [a for a in list(self._axes.lines) + list(self._axes.texts) + if a not in before] + + # Axes.add_artist clips to the axes patch but the exterior + # labels sit just outside it. + for artist in self.artists: + if isinstance(artist, matplotlib.text.Text): + artist.set_clip_on(False) + + def _removeArtists(self): + """Remove the AST artists from the axes.""" + for artist in self.artists: + try: + artist.remove() + except (ValueError, NotImplementedError): + # Already gone, e.g. the axes were cleared. + pass + self.artists = [] + + def remove(self): + """Remove the AST axes and restore matplotlib's axis furniture.""" + self._removeArtists() + + xVisible, yVisible, spines = self._savedVisibility + self._axes.xaxis.set_visible(xVisible) + self._axes.yaxis.set_visible(yVisible) + for name, visible in spines.items(): + self._axes.spines[name].set_visible(visible) +``` + +- [ ] **Step 4: Run tests to verify they pass** + +Same command as Step 2, plus the Task 1 test: + +```bash +pytest tests/test_display_matplotlib.py -v +``` + +Expected: all PASS. + +- [ ] **Step 5: Run flake8 and commit** + +```bash +flake8 python tests +git add python/lsst/display/matplotlib/wcsAxes.py tests/test_display_matplotlib.py +git commit -m "Add WcsAxesManager to draw AST sky axes on matplotlib axes + +Co-Authored-By: Claude Fable 5 " +``` + +--- + +### Task 3: Redraw on view change and sexagesimal toggle + +**Files:** +- Modify: `python/lsst/display/matplotlib/wcsAxes.py` +- Test: `tests/test_display_matplotlib.py` + +**Interfaces:** +- Consumes: `WcsAxesManager` from Task 2. +- Produces: `WcsAxesManager.setSexagesimal(useSexagesimal)`; automatic redraw when the axes' x or y limits change; `remove()` additionally disconnects the callbacks. Task 4/5 rely on both. + +- [ ] **Step 1: Write the failing tests** + +Add to `WcsAxesManagerTestCase` in `tests/test_display_matplotlib.py`: + +```python + def testRedrawOnLimitChange(self): + ax, manager = makeWcsAxes() + oldArtists = list(manager.artists) + ax.set_xlim(50, 150) + ax.set_ylim(40, 110) + self.assertGreater(len(manager.artists), 0) + self.assertNotEqual(manager.artists, oldArtists) + # The old artists must be gone from the axes. + axesArtists = set(ax.lines) | set(ax.texts) + for artist in oldArtists: + self.assertNotIn(artist, axesArtists) + + def testNoRedrawAfterRemove(self): + ax, manager = makeWcsAxes() + manager.remove() + ax.set_xlim(50, 150) + self.assertEqual(manager.artists, []) + self.assertEqual(len(ax.lines), 0) + + def testSetSexagesimal(self): + ax, manager = makeWcsAxes() + + def tickLabels(): + return [a.get_text() for a in manager.artists + if isinstance(a, matplotlib.text.Text)] + + self.assertFalse(any(":" in label for label in tickLabels())) + manager.setSexagesimal(True) + self.assertTrue(any(":" in label for label in tickLabels())) + manager.setSexagesimal(False) + self.assertFalse(any(":" in label for label in tickLabels())) +``` + +- [ ] **Step 2: Run tests to verify they fail** + +```bash +export LSSTSW=/Users/timj/work/lsstsw +source "$LSSTSW/bin/envconfig" +setup -t b8267 lsst_distrib +cd /Users/timj/work/lsstsw/build/display_matplotlib +setup -k -r . +pytest tests/test_display_matplotlib.py::WcsAxesManagerTestCase -v +``` + +Expected: the three new tests FAIL (`testRedrawOnLimitChange` on `assertNotEqual`, `testSetSexagesimal` with `AttributeError: ... has no attribute 'setSexagesimal'`). + +- [ ] **Step 3: Implement** + +In `WcsAxesManager.__init__`, after the spine-hiding loop and before `self.draw()`, add: + +```python + self._redrawing = False + self._callbackIds = [ + axes.callbacks.connect("xlim_changed", self._onLimitsChanged), + axes.callbacks.connect("ylim_changed", self._onLimitsChanged), + ] +``` + +Add these methods to `WcsAxesManager`: + +```python + def _onLimitsChanged(self, axes): + """Redraw for new view limits; matplotlib callback.""" + if self._redrawing: + return + self._redrawing = True + try: + self.draw() + finally: + self._redrawing = False + + def setSexagesimal(self, useSexagesimal): + """Switch between sexagesimal and decimal degree labels. + + Parameters + ---------- + useSexagesimal : `bool` + Format labels as e.g. HH:MM:SS.s iff True. + """ + useSexagesimal = bool(useSexagesimal) + if useSexagesimal != self._useSexagesimal: + self._useSexagesimal = useSexagesimal + self.draw() +``` + +In `remove()`, before `self._removeArtists()`, add: + +```python + for callbackId in self._callbackIds: + self._axes.callbacks.disconnect(callbackId) + self._callbackIds = [] +``` + +- [ ] **Step 4: Run tests to verify they pass** + +```bash +pytest tests/test_display_matplotlib.py -v +``` + +Expected: all PASS. + +- [ ] **Step 5: Run flake8 and commit** + +```bash +flake8 python tests +git add python/lsst/display/matplotlib/wcsAxes.py tests/test_display_matplotlib.py +git commit -m "Redraw WCS axes on view changes and support sexagesimal toggle + +Co-Authored-By: Claude Fable 5 " +``` + +--- + +### Task 4: DisplayImpl integration + +**Files:** +- Modify: `python/lsst/display/matplotlib/matplotlib.py` +- Test: `tests/test_display_matplotlib.py` + +**Interfaces:** +- Consumes: `astFrameSetFromWcs`, `WcsAxesManager` (Tasks 1-3); `DisplayImpl.__init__`, `_mtv` (existing). +- Produces: `DisplayImpl` kwargs `useWcsAxes=True`, `wcsGrid=False`, `astPlotOptions=""` reachable via `afwDisplay.Display(...)`; `self._wcsAxesManagers` dict mapping `matplotlib.axes.Axes` to `WcsAxesManager`; method `_drawWcsAxes(ax, wcs)`. Task 5 uses `self._wcsAxesManagers`. + +- [ ] **Step 1: Write the failing tests** + +In `tests/test_display_matplotlib.py`, add to the import block (with the other `lsst` imports): + +```python +import lsst.afw.image as afwImage # noqa: E402 +``` + +Then add: + +```python +class WcsAxesDisplayTestCase(lsst.utils.tests.TestCase): + """Tests for WCS axes at the afwDisplay level.""" + + frame = 100 # keep test figures separate from any other test + + def setUp(self): + afwDisplay.setDefaultBackend("matplotlib") + # A distinct frame per test so each test gets a fresh figure. + WcsAxesDisplayTestCase.frame += 1 + self.frame = WcsAxesDisplayTestCase.frame + + def tearDown(self): + afwDisplay.delAllDisplays() + + @staticmethod + def makeExposure(): + exposure = afwImage.ExposureF(200, 150) + exposure.image.array[:] = 1.0 + exposure.setWcs(makeTestWcs()) + return exposure + + def testWcsAxesByDefault(self): + display = afwDisplay.Display(frame=self.frame) + display.mtv(self.makeExposure()) + impl = display._impl + ax = impl._figure.gca() + self.assertIn(ax, impl._wcsAxesManagers) + self.assertGreater(len(impl._wcsAxesManagers[ax].artists), 0) + self.assertFalse(ax.xaxis.get_visible()) + + def testUseWcsAxesFalse(self): + display = afwDisplay.Display(frame=self.frame, useWcsAxes=False) + display.mtv(self.makeExposure()) + impl = display._impl + ax = impl._figure.gca() + self.assertEqual(impl._wcsAxesManagers, {}) + self.assertTrue(ax.xaxis.get_visible()) + + def testNoWcsMeansPixelAxes(self): + display = afwDisplay.Display(frame=self.frame) + display.mtv(afwImage.ImageF(200, 150)) + impl = display._impl + ax = impl._figure.gca() + self.assertEqual(impl._wcsAxesManagers, {}) + self.assertTrue(ax.xaxis.get_visible()) + + def testWcsGridOption(self): + display = afwDisplay.Display(frame=self.frame, wcsGrid=True) + display.mtv(self.makeExposure()) + impl = display._impl + ax = impl._figure.gca() + manager = impl._wcsAxesManagers[ax] + self.assertIn("Grid=1", manager._plotOptions()) + + def testZoomRedraws(self): + display = afwDisplay.Display(frame=self.frame) + display.mtv(self.makeExposure()) + impl = display._impl + ax = impl._figure.gca() + manager = impl._wcsAxesManagers[ax] + oldArtists = list(manager.artists) + display.zoom(4) + display.pan(100, 75) + self.assertGreater(len(manager.artists), 0) + self.assertNotEqual(manager.artists, oldArtists) +``` + +- [ ] **Step 2: Run tests to verify they fail** + +```bash +export LSSTSW=/Users/timj/work/lsstsw +source "$LSSTSW/bin/envconfig" +setup -t b8267 lsst_distrib +cd /Users/timj/work/lsstsw/build/display_matplotlib +setup -k -r . +pytest tests/test_display_matplotlib.py::WcsAxesDisplayTestCase -v +``` + +Expected: FAIL with `AttributeError: 'DisplayImpl' object has no attribute '_wcsAxesManagers'` (and `TypeError` for the kwarg tests). + +- [ ] **Step 3: Implement** + +In `python/lsst/display/matplotlib/matplotlib.py`: + +**(a)** Add to the import block (after `import lsst.geom as geom`, around line 51): + +```python +from .wcsAxes import WcsAxesManager, astFrameSetFromWcs +``` + +**(b)** Change the `__init__` signature (line 100-102) to add the three kwargs: + +```python + def __init__(self, display, verbose=False, + interpretMaskBits=True, mtvOrigin=afwImage.PARENT, fastMaskDisplay=False, + reopenPlot=False, useSexagesimal=False, dpi=None, + useWcsAxes=True, wcsGrid=False, astPlotOptions="", *args, **kwargs): +``` + +**(c)** Extend the `__init__` docstring parameter list (after the `@param dpi` line): + +``` + @param useWcsAxes If True (the default) and the data + have a WCS, draw sky coordinate axes + using AST instead of pixel axes + @param wcsGrid If True draw the full curvilinear + coordinate grid across the image + (only used with useWcsAxes) + @param astPlotOptions Extra AST Plot attribute settings, + e.g. "Colour(grid)=red"; these + override the defaults +``` + +**(d)** In the `__init__` body, after `self._mtvOrigin = mtvOrigin` (line 158), add: + +```python + self._useWcsAxes = useWcsAxes + self._wcsGrid = wcsGrid + self._astPlotOptions = astPlotOptions + self._wcsAxesManagers = {} # matplotlib Axes -> WcsAxesManager +``` + +**(e)** In `_mtv`, replace (lines 336-337): + +```python + ax = self._figure.gca() + ax.cla() +``` + +with: + +```python + ax = self._figure.gca() + manager = self._wcsAxesManagers.pop(ax, None) + if manager is not None: + manager.remove() + ax.cla() +``` + +**(f)** In `_mtv`, after the `self._title = title` line (line 349), add: + +```python + if wcs is not None and self._useWcsAxes: + self._drawWcsAxes(ax, wcs) +``` + +**(g)** Add this method after `_i_mtv` (after line 491): + +```python + def _drawWcsAxes(self, ax, wcs): + """Draw sky coordinate axes with AST, hiding the pixel axes. + + Falls back to ordinary pixel axes with a warning if AST cannot + draw the given WCS. + """ + # Freeze the view at the image extent; autoscaling in response + # to AST artists drawn outside the view would change the limits + # and retrigger drawing. + ax.set_xlim(*ax.get_xlim()) + ax.set_ylim(*ax.get_ylim()) + + try: + frameSet = astFrameSetFromWcs(wcs) + self._wcsAxesManagers[ax] = WcsAxesManager( + ax, frameSet, wcs.getPixelScale().asDegrees(), + grid=self._wcsGrid, + useSexagesimal=self._useSexagesimal[0], + extraOptions=self._astPlotOptions) + except Exception as e: + print(f"Unable to draw WCS axes ({e}); using pixel coordinates", + file=sys.stderr) +``` + +- [ ] **Step 4: Run tests to verify they pass** + +```bash +pytest tests/test_display_matplotlib.py -v +``` + +Expected: all PASS. + +- [ ] **Step 5: Run flake8 and commit** + +```bash +flake8 python tests +git add python/lsst/display/matplotlib/matplotlib.py tests/test_display_matplotlib.py +git commit -m "Draw WCS axes by default when displaying data with a WCS + +Co-Authored-By: Claude Fable 5 " +``` + +--- + +### Task 5: Lifecycle — erase, close, useSexagesimal + +**Files:** +- Modify: `python/lsst/display/matplotlib/matplotlib.py` +- Test: `tests/test_display_matplotlib.py` + +**Interfaces:** +- Consumes: `self._wcsAxesManagers` and `WcsAxesManager.artists`/`remove()`/`setSexagesimal()` (Tasks 2-4). +- Produces: `_erase()` that clears overlay glyphs but keeps WCS axes (and works on matplotlib ≥ 3.7, where the current implementation raises `AttributeError`); `_close()` that releases managers; `useSexagesimal()` that updates the axis labels live. + +- [ ] **Step 1: Write the failing tests** + +Add to `WcsAxesDisplayTestCase`: + +```python + def testErasePreservesWcsAxes(self): + display = afwDisplay.Display(frame=self.frame) + display.mtv(self.makeExposure()) + impl = display._impl + ax = impl._figure.gca() + manager = impl._wcsAxesManagers[ax] + nAstLines = sum(1 for a in manager.artists if a in set(ax.lines)) + display.dot("+", 100, 75) + self.assertGreater(len(ax.lines), nAstLines) + display.erase() + self.assertEqual(sum(1 for a in ax.lines), nAstLines) + axesArtists = set(ax.lines) | set(ax.texts) + for artist in manager.artists: + self.assertIn(artist, axesArtists) + + def testRepeatedMtvDoesNotLeak(self): + display = afwDisplay.Display(frame=self.frame) + display.mtv(self.makeExposure()) + impl = display._impl + ax = impl._figure.gca() + nTexts = len(ax.texts) + display.mtv(self.makeExposure()) + ax = impl._figure.gca() + self.assertEqual(len(impl._wcsAxesManagers), 1) + self.assertEqual(len(ax.texts), nTexts) + + def testUseSexagesimalRedraws(self): + display = afwDisplay.Display(frame=self.frame) + display.mtv(self.makeExposure()) + impl = display._impl + ax = impl._figure.gca() + manager = impl._wcsAxesManagers[ax] + + def hasColon(): + return any(":" in a.get_text() for a in manager.artists + if isinstance(a, matplotlib.text.Text)) + + self.assertFalse(hasColon()) + display.useSexagesimal(True) + self.assertTrue(hasColon()) + + def testCloseReleasesManagers(self): + display = afwDisplay.Display(frame=self.frame) + display.mtv(self.makeExposure()) + impl = display._impl + impl._close() + self.assertEqual(impl._wcsAxesManagers, {}) +``` + +Note: `display.useSexagesimal(True)` reaches the impl through +`lsst.afw.display.interface.Display.__getattr__` (interface.py:222), +which delegates unknown attributes to the impl, the same way +`display.savefig` works today (verified in the installed afw). + +- [ ] **Step 2: Run tests to verify they fail** + +```bash +export LSSTSW=/Users/timj/work/lsstsw +source "$LSSTSW/bin/envconfig" +setup -t b8267 lsst_distrib +cd /Users/timj/work/lsstsw/build/display_matplotlib +setup -k -r . +pytest tests/test_display_matplotlib.py::WcsAxesDisplayTestCase -v +``` + +Expected: `testErasePreservesWcsAxes` FAILS with `AttributeError: property 'lines' of 'Axes' object has no setter` (the pre-existing `_erase` bug); `testUseSexagesimalRedraws` fails on the label assertion; `testCloseReleasesManagers` may fail on leftover managers. + +- [ ] **Step 3: Implement** + +In `python/lsst/display/matplotlib/matplotlib.py`: + +**(a)** Replace `_erase` (lines 537-544): + +```python + def _erase(self): + """Erase the display, keeping any WCS axes""" + + for axis in self._figure.axes: + manager = self._wcsAxesManagers.get(axis) + keep = set(manager.artists) if manager is not None else set() + for artist in list(axis.lines) + list(axis.texts): + if artist not in keep: + artist.remove() + + self._figure.canvas.draw_idle() +``` + +**(b)** In `_close` (lines 176-181), add before `self._image = None`: + +```python + for manager in self._wcsAxesManagers.values(): + manager.remove() + self._wcsAxesManagers = {} +``` + +**(c)** In `useSexagesimal` (line 268), after `self._useSexagesimal[0] = useSexagesimal`, add: + +```python + for manager in self._wcsAxesManagers.values(): + manager.setSexagesimal(useSexagesimal) + self._figure.canvas.draw_idle() +``` + +- [ ] **Step 4: Run tests to verify they pass** + +```bash +pytest tests/test_display_matplotlib.py -v +``` + +Expected: all PASS. + +- [ ] **Step 5: Run flake8 and commit** + +```bash +flake8 python tests +git add python/lsst/display/matplotlib/matplotlib.py tests/test_display_matplotlib.py +git commit -m "Preserve WCS axes across erase and update them on option changes + +Co-Authored-By: Claude Fable 5 " +``` + +--- + +### Task 6: Full-suite verification and visual check + +**Files:** +- No production changes expected; fix anything the checks surface. + +- [ ] **Step 1: Run the complete test suite and flake8** + +```bash +export LSSTSW=/Users/timj/work/lsstsw +source "$LSSTSW/bin/envconfig" +setup -t b8267 lsst_distrib +cd /Users/timj/work/lsstsw/build/display_matplotlib +setup -k -r . +pytest tests/test_display_matplotlib.py -v +flake8 python tests +``` + +Expected: all tests PASS; flake8 silent. + +- [ ] **Step 2: Visual verification** + +Write this script to the session scratchpad directory as `visual_check.py` and run it there (same EUPS setup block first): + +```python +import matplotlib +matplotlib.use("Agg") + +import lsst.afw.display as afwDisplay +import lsst.afw.geom as afwGeom +import lsst.afw.image as afwImage +import lsst.geom +import numpy as np + +afwDisplay.setDefaultBackend("matplotlib") + +exposure = afwImage.ExposureF(200, 150) +exposure.image.array[:] = np.random.default_rng(42).normal( + 100, 10, exposure.image.array.shape) +exposure.setWcs(afwGeom.makeSkyWcs( + crpix=lsst.geom.Point2D(100, 75), + crval=lsst.geom.SpherePoint(30.0, -45.0, lsst.geom.degrees), + cdMatrix=afwGeom.makeCdMatrix(scale=3.0*lsst.geom.arcseconds, + orientation=30*lsst.geom.degrees))) + +cases = [ + ("check_default.png", dict()), + ("check_grid.png", dict(wcsGrid=True)), + ("check_sexagesimal.png", dict(useSexagesimal=True)), + ("check_pixel.png", dict(useWcsAxes=False)), +] +for frame, (filename, kwargs) in enumerate(cases, start=1): + display = afwDisplay.Display(frame=frame, **kwargs) + display.scale("linear", "minmax") + display.mtv(exposure, title=filename) + display._impl.savefig(filename, dpi=100) + +display = afwDisplay.Display(frame=10) +display.mtv(exposure, title="zoomed") +display.zoom(8) +display.pan(100, 75) +display._impl.savefig("check_zoom.png", dpi=100) +print("wrote all PNGs") +``` + +- [ ] **Step 3: Inspect each PNG with the Read tool and confirm** + +- `check_default.png`: sky tick labels in decimal degrees on all edges, "Right ascension"/"Declination" labels, no matplotlib pixel ticks, title shown. +- `check_grid.png`: curvilinear grid lines crossing the image at ~30° (the WCS rotation). +- `check_sexagesimal.png`: labels like `2:00:00` / `-45:00`. +- `check_pixel.png`: plain matplotlib pixel axes (numbers 0-200), no sky labels. +- `check_zoom.png`: labels present and consistent with a ~25-pixel-wide view (more decimal places than the default view). + +- [ ] **Step 4: Commit any fixes** + +If steps 1-3 surfaced fixes, run flake8 and commit them: + +```bash +flake8 python tests +git add -A python tests +git commit -m "Fix issues found during WCS axes verification + +Co-Authored-By: Claude Fable 5 " +``` diff --git a/docs/superpowers/specs/2026-07-03-wcs-axes-design.md b/docs/superpowers/specs/2026-07-03-wcs-axes-design.md new file mode 100644 index 0000000..76c2f29 --- /dev/null +++ b/docs/superpowers/specs/2026-07-03-wcs-axes-design.md @@ -0,0 +1,99 @@ +# Design: WCS axes for display_matplotlib via starlink-pyast + +## Goal + +When `mtv()`/`image()` displays data that has a WCS, the axes show sky coordinates (RA/Dec) instead of pixel coordinates. +All axis furniture — border, ticks, tick labels, axis labels, and optional grid — is drawn by AST through pyast's matplotlib grf plugin. +The matplotlib `Axes` object remains as the container: it hosts `imshow` and overlay artists, and its data coordinates are still LSST pixel coordinates. +Its own spines, ticks, and labels are switched off entirely; AST is the only thing drawing visible axes. + +This uses the AST-native representation of the WCS throughout, with no FITS approximation. + +## Usage + +```python +afwDisplay.setDefaultBackend("matplotlib") +display = afwDisplay.Display() # WCS axes on by default when a WCS is present +display.scale('asinh', 'zscale') +display.image(cutout) +pyplot.show() +``` + +`afwDisplay.Display(frame=None, backend=None, **kwargs)` forwards `**kwargs` to the backend `DisplayImpl` constructor, which is how the new options reach this plugin (the same mechanism as the existing `useSexagesimal`, `fastMaskDisplay`, and `dpi` options). + +## Architecture and data flow + +### WCS conversion + +A module-level helper converts an afw `SkyWcs` to a pyast `FrameSet`: + +1. Write `wcs.getFrameDict()` through an `astshim.Channel` into an `astshim.StringStream`. +2. Read the resulting AST native text back with a `starlink.Ast.Channel`, using a small adapter class that exposes an `astsource()` method returning one line per call. + +This is AST-native serialization, so the round-trip is exact. +This was verified in the b8267 stack environment: the resulting FrameSet's `tran()` agrees with `wcs.pixelToSky()` to the last bit. +The FrameSet arrives with base frame `PIXELS` in LSST 0-based pixel coordinates — identical to the axes data coordinates established by `imshow(extent=...)` — and current frame `SKY`. +No coordinate offsets are needed. + +### Drawing + +In `_mtv`, when a WCS is present and WCS axes are enabled, after the image (and mask) are drawn: + +1. Switch off all matplotlib axis furniture on the axes: spines, ticks, tick labels, and axis labels. +2. Create `Ast.Plot(frameset, gbox, bbox, options, grf_matplotlib(ax))`, where `gbox` and `bbox` are both the current view limits in pixel coordinates (they coincide because the data coordinates are pixel coordinates). +3. Call `plot.grid()`; AST draws the border, ticks, tick labels, axis labels, and (optionally) the curvilinear grid as matplotlib artists on the axes. +4. Track the artists AST created so they can be removed cleanly later. + +### Lifecycle + +The AST axes are rebuilt — old artists removed, a new `Ast.Plot` created for the new view — whenever the view changes: + +- `xlim_changed`/`ylim_changed` axes callbacks cover both interactive toolbar navigation and programmatic `display.zoom()`/`display.pan()` (which set axis limits). + The two callbacks are coalesced and guarded against re-entry so one zoom triggers one rebuild. +- `_erase()` must not destroy the WCS axes. + Its current behavior of clearing all lines and texts would delete the AST artists, so erase instead removes only artists that are not in the tracked AST set. +- A subsequent `mtv()` rebuilds from scratch without leaking artists or callbacks. +- `_close()` drops the FrameSet and disconnects the callbacks. + +When there is no WCS, or WCS axes are disabled, behavior is exactly the current pixel-axes behavior and nothing is hidden. + +## API and formatting + +New constructor kwargs on the backend (reachable via `afwDisplay.Display(...)`): + +- `useWcsAxes=True` — draw AST sky axes when the displayed data has a WCS; `False` restores plain pixel axes. +- `wcsGrid=False` — when `True`, draw the full curvilinear grid across the image (AST `Grid=1`); the default is edge ticks and labels only. +- `astPlotOptions=""` — an AST attribute string (e.g. `"Colour(grid)=red, Size(TextLab)=12"`) appended to the Plot options the backend constructs. + Because it comes last, it overrides the defaults; this is the escape hatch for styling without growing a zoo of kwargs. + +Label formatting follows the existing `useSexagesimal` flag: + +- `False` (default): decimal degrees, set via the SkyFrame `Format` attributes on the Plot. +- `True`: AST's native hms/dms sky formatting. +- Toggling `display.useSexagesimal(...)` at runtime rebuilds the AST axes so the axis labels stay consistent with the cursor readout. + +Axis label text comes from the SkyFrame (`Label(1)`/`Label(2)`, e.g. "Right ascension" / "Declination"), not hardcoded strings. + +## Edge cases + +- User-supplied `Figure` objects and the `fig.sca(axis)` multi-panel pattern keep working; AST state (artists, callbacks, FrameSet) is tracked per-axes. +- Mask-only and WCS-less displays follow the unchanged pixel-axes path. +- The colorbar (`make_axes_locatable`) and matplotlib's `set_title` are unaffected; the title remains matplotlib's. +- `starlink.Ast` and `starlink.Grf` are imported at module top level; starlink-pyast becomes a hard dependency of display_matplotlib (it is provided by rubin-env). + +## Error handling + +If `Ast.Plot` construction or `plot.grid()` raises for a pathological WCS, emit a warning and fall back to pixel axes rather than failing `mtv`. + +## Testing + +Extend `tests/test_display_matplotlib.py` (Agg backend): + +- `mtv` of an exposure with a WCS produces AST artists and hides the matplotlib ticks and spines; axis labels are present. +- The SkyWcs → pyast FrameSet round-trip helper matches `pixelToSky` for sample points. +- `useWcsAxes=False` reproduces the existing pixel-axes behavior. +- `zoom()` rebuilds the AST artists for the new view. +- `erase()` preserves the WCS axes while clearing overlay glyphs. +- Repeated `mtv()` calls do not leak artists or callbacks. + +Manual verification in the EUPS stack environment with saved PNGs: decimal and sexagesimal labels, grid on and off, zoomed and unzoomed views. diff --git a/python/lsst/display/matplotlib/matplotlib.py b/python/lsst/display/matplotlib/matplotlib.py index 58d51ca..319030b 100644 --- a/python/lsst/display/matplotlib/matplotlib.py +++ b/python/lsst/display/matplotlib/matplotlib.py @@ -50,6 +50,8 @@ import lsst.afw.geom as afwGeom import lsst.geom as geom +from .wcsAxes import WcsAxesManager, astFrameSetFromWcs + # # Set the list of backends which support _getEvent and thus interact() # @@ -99,7 +101,8 @@ class DisplayImpl(virtualDevice.DisplayImpl): """ def __init__(self, display, verbose=False, interpretMaskBits=True, mtvOrigin=afwImage.PARENT, fastMaskDisplay=False, - reopenPlot=False, useSexagesimal=False, dpi=None, *args, **kwargs): + reopenPlot=False, useSexagesimal=False, dpi=None, + useWcsAxes=True, wcsGrid=False, astPlotOptions="", *args, **kwargs): """ Initialise a matplotlib display @@ -118,6 +121,19 @@ def __init__(self, display, verbose=False, May be changed by calling display.useSexagesimal() @param dpi Number of dpi (passed to pyplot.figure) + @param useWcsAxes If True (the default) and the data + have a WCS, draw sky coordinate axes + using AST instead of pixel axes + @param wcsGrid If True draw the full curvilinear + coordinate grid across the image + (only used with useWcsAxes) + @param astPlotOptions Extra AST Plot attribute settings, + e.g. "Colour(grid)=2" for a red + grid; these override the defaults. + N.b. AST colours are 1-based grf + colour table indices, not names + (see lsst.display.matplotlib + .wcsAxes.WcsAxesManager) The `frame` argument to `Display` may be a matplotlib figure; this permits code such as @@ -156,6 +172,10 @@ def __init__(self, display, verbose=False, self._fastMaskDisplay = fastMaskDisplay self._useSexagesimal = [useSexagesimal] # use an array so we can modify the value in format_coord self._mtvOrigin = mtvOrigin + self._useWcsAxes = useWcsAxes + self._wcsGrid = wcsGrid + self._astPlotOptions = astPlotOptions + self._wcsAxesManagers = {} # matplotlib Axes -> WcsAxesManager self._mappable_ax = None self._colorbar_ax = None self._image_colormap = matplotlib.cm.gray @@ -175,6 +195,9 @@ def __init__(self, display, verbose=False, def _close(self): """!Close the display, cleaning up any allocated resources""" + for manager in self._wcsAxesManagers.values(): + manager.remove() + self._wcsAxesManagers = {} self._image = None self._mask = None self._wcs = None @@ -279,6 +302,10 @@ def useSexagesimal(self, useSexagesimal): """Are we formatting coordinates as HH:MM:SS.ss?""" self._useSexagesimal[0] = useSexagesimal + for manager in self._wcsAxesManagers.values(): + manager.setSexagesimal(useSexagesimal) + self._figure.canvas.draw_idle() + def wait(self, prompt="[c(ontinue) p(db)] :", allowPdb=True): """Wait for keyboard input @@ -334,6 +361,9 @@ def _mtv(self, image, mask=None, wcs=None, title="", metadata=None): self._scaleArgs['unit'], *self._scaleArgs['args'], **self._scaleArgs['kwargs']) ax = self._figure.gca() + manager = self._wcsAxesManagers.pop(ax, None) + if manager is not None: + manager.remove() ax.cla() self._i_mtv(image, wcs, title, False) @@ -348,6 +378,9 @@ def _mtv(self, image, mask=None, wcs=None, title="", metadata=None): self._title = title + if wcs is not None and self._useWcsAxes: + self._drawWcsAxes(ax, wcs) + def format_coord(x, y, wcs=self._wcs, x0=self._xy0[0], y0=self._xy0[1], origin=afwImage.PARENT, bbox=self._image.getBBox(afwImage.PARENT), _useSexagesimal=self._useSexagesimal): @@ -490,6 +523,29 @@ def _i_mtv(self, data, wcs, title, isMask): self._figure.canvas.draw_idle() + def _drawWcsAxes(self, ax, wcs): + """Draw sky coordinate axes with AST, hiding the pixel axes. + + Falls back to ordinary pixel axes with a warning if AST cannot + draw the given WCS. + """ + # Freeze the view at the image extent; autoscaling in response + # to AST artists drawn outside the view would change the limits + # and retrigger drawing. + ax.set_xlim(*ax.get_xlim()) + ax.set_ylim(*ax.get_ylim()) + + try: + frameSet = astFrameSetFromWcs(wcs) + self._wcsAxesManagers[ax] = WcsAxesManager( + ax, frameSet, + grid=self._wcsGrid, + useSexagesimal=self._useSexagesimal[0], + extraOptions=self._astPlotOptions) + except Exception as e: + print(f"Unable to draw WCS axes ({e}); using pixel coordinates", + file=sys.stderr) + def _i_setImage(self, image, mask=None, wcs=None): """Save the current image, mask, wcs, and XY0""" self._image = image @@ -535,11 +591,14 @@ def _flush(self): pass def _erase(self): - """Erase the display""" + """Erase the display, keeping any WCS axes""" for axis in self._figure.axes: - axis.lines = [] - axis.texts = [] + manager = self._wcsAxesManagers.get(axis) + keep = set(manager.artists) if manager is not None else set() + for artist in list(axis.lines) + list(axis.texts) + list(axis.patches): + if artist not in keep: + artist.remove() self._figure.canvas.draw_idle() diff --git a/python/lsst/display/matplotlib/wcsAxes.py b/python/lsst/display/matplotlib/wcsAxes.py new file mode 100644 index 0000000..2c77f0c --- /dev/null +++ b/python/lsst/display/matplotlib/wcsAxes.py @@ -0,0 +1,281 @@ +# +# LSST Data Management System +# Copyright 2026 LSST Corporation. +# +# This product includes software developed by the +# LSST Project (http://www.lsst.org/). +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the LSST License Statement and +# the GNU General Public License along with this program. If not, +# see . +# + +"""Draw sky coordinate axes on a matplotlib Axes using AST. + +The axes are drawn by `starlink.Ast.Plot` using the matplotlib grf +plugin, so arbitrary AST WCS transforms are supported with no FITS +approximation. +""" + +__all__ = ["astFrameSetFromWcs", "WcsAxesManager"] + +import astshim +import matplotlib.text +import starlink.Ast +from matplotlib.backend_bases import TimerBase +from starlink.Grf import grf_matplotlib + + +class _AstStringSource: + """Present AST native serialization text to `starlink.Ast.Channel`. + + `starlink.Ast.Channel` reads via an object with an ``astsource`` + method returning one line per call, and `None` at end of input. + + Parameters + ---------- + text : `str` + AST native serialization of an AST object. + """ + + def __init__(self, text): + self._lines = text.splitlines() + self._pos = 0 + + def astsource(self): + if self._pos >= len(self._lines): + return None + line = self._lines[self._pos] + self._pos += 1 + return line + + +def astFrameSetFromWcs(wcs): + """Convert an afw SkyWcs to a starlink-pyast FrameSet. + + Parameters + ---------- + wcs : `lsst.afw.geom.SkyWcs` + The WCS to convert. + + Returns + ------- + frameSet : `starlink.Ast.FrameSet` + FrameSet whose base frame is LSST 0-based pixel coordinates + (domain ``PIXELS``) and whose current frame is sky coordinates + in radians (domain ``SKY``). + + Notes + ----- + The conversion serializes the WCS's underlying `astshim.FrameDict` + to AST native text and reads it back with pyast, so the result is + exact; no FITS approximation is involved. + """ + stream = astshim.StringStream() + astshim.Channel(stream).write(wcs.getFrameDict()) + channel = starlink.Ast.Channel(_AstStringSource(stream.getSinkData())) + return channel.read() + + +class WcsAxesManager: + """Manage AST-drawn sky coordinate axes on a matplotlib Axes. + + All visible axis furniture (border, ticks, labels, optional grid) + is drawn by AST; matplotlib's own axis furniture is hidden while + this manager is active and restored by `remove`. + + The axes' data coordinates must be LSST 0-based pixel coordinates + matching the base frame of ``frameSet``, as produced by + `astFrameSetFromWcs`. + + Parameters + ---------- + axes : `matplotlib.axes.Axes` + The axes to draw on. + frameSet : `starlink.Ast.FrameSet` + Pixel-to-sky transform, base frame in axes data coordinates. + grid : `bool`, optional + Draw the full curvilinear coordinate grid across the image? + useSexagesimal : `bool`, optional + Format labels as e.g. HH:MM:SS.s rather than decimal degrees. + extraOptions : `str`, optional + Comma-separated AST Plot attribute settings appended after the + defaults, so they take precedence (e.g. + ``"Colour(grid)=2, Width(border)=2"``). + + Notes + ----- + AST colour values are 1-based indices into the grf colour table + (see `starlink.Grf.grf_matplotlib`), not colour names: + 1=default, 2=red, 3=green, 4=blue, 5=cyan, 6=magenta, 7=yellow, + 8=black, 9=dark grey, 10=grey, 11=light grey, 12=white. + """ + + def __init__(self, axes, frameSet, *, + grid=False, useSexagesimal=False, extraOptions=""): + self._axes = axes + self._frameSet = frameSet + self._grid = grid + self._useSexagesimal = useSexagesimal + self._extraOptions = extraOptions + self.artists = [] + + self._savedVisibility = ( + axes.xaxis.get_visible(), + axes.yaxis.get_visible(), + {name: spine.get_visible() for name, spine in axes.spines.items()}, + ) + axes.xaxis.set_visible(False) + axes.yaxis.set_visible(False) + for spine in axes.spines.values(): + spine.set_visible(False) + + self._redrawing = False + self._redrawTimer = None + self._callbackIds = [ + axes.callbacks.connect("xlim_changed", self._onLimitsChanged), + axes.callbacks.connect("ylim_changed", self._onLimitsChanged), + ] + + try: + self.draw() + except Exception: + # Restore the matplotlib furniture so the axes remain + # usable as ordinary pixel axes. + self.remove() + raise + + def _plotOptions(self): + """Return the AST Plot options string for the current state.""" + options = ["DrawTitle=0"] + options.append("Grid=1" if self._grid else "Grid=0") + if not self._useSexagesimal: + # Decimal degrees; ".*" lets the Plot choose the precision + # needed to keep adjacent labels distinct. + options.append("Format(1)=d.*") + options.append("Format(2)=d.*") + if self._extraOptions: + options.append(self._extraOptions) + return ", ".join(options) + + def draw(self): + """Draw (or redraw) the AST axes for the current view limits.""" + self._removeArtists() + + xlo, xhi = self._axes.get_xlim() + ylo, yhi = self._axes.get_ylim() + box = (xlo, ylo, xhi, yhi) + + before = set(self._axes.lines) | set(self._axes.texts) + try: + plot = starlink.Ast.Plot(self._frameSet, box, box, + grf_matplotlib(self._axes), + self._plotOptions()) + plot.grid() + except Exception: + # grid() may have drawn some artists before failing. + for artist in list(self._axes.lines) + list(self._axes.texts): + if artist not in before: + artist.remove() + raise + self.artists = [a for a in list(self._axes.lines) + list(self._axes.texts) + if a not in before] + + # Axes.add_artist clips to the axes patch but the exterior + # labels sit just outside it. + for artist in self.artists: + if isinstance(artist, matplotlib.text.Text): + artist.set_clip_on(False) + + def _removeArtists(self): + """Remove the AST artists from the axes.""" + for artist in self.artists: + try: + artist.remove() + except (ValueError, NotImplementedError): + # Already gone, e.g. the axes were cleared. + pass + self.artists = [] + + # Idle time before redrawing after a view change. Interactive + # panning and zooming change the limits on every mouse-motion + # event, and a full AST rebuild per event is far too slow. + _redrawDelayMs = 200 + + def _onLimitsChanged(self, axes): + """Schedule a redraw for new view limits; matplotlib callback. + + The redraw is debounced with a one-shot timer so that a stream + of limit changes (interactive panning or zooming, or the x/y + pair from a single zoom) produces a single rebuild once the + view settles. Backends without a running event loop cannot + fire timers, so there the redraw happens immediately. + """ + if self._redrawing: + return + if self._redrawTimer is None: + canvas = self._axes.get_figure().canvas + timer = canvas.new_timer(interval=self._redrawDelayMs) + if type(timer) is TimerBase: + # The base class is a no-op that never fires. + self._redrawNow() + return + timer.single_shot = True + timer.add_callback(self._redrawNow) + self._redrawTimer = timer + self._redrawTimer.stop() + self._redrawTimer.start() + + def _redrawNow(self): + """Rebuild the AST axes for the current view and repaint.""" + if not self._callbackIds: + # The manager was removed while a redraw was pending. + return + if self._redrawing: + return + self._redrawing = True + try: + self.draw() + finally: + self._redrawing = False + self._axes.get_figure().canvas.draw_idle() + + def setSexagesimal(self, useSexagesimal): + """Switch between sexagesimal and decimal degree labels. + + Parameters + ---------- + useSexagesimal : `bool` + Format labels as e.g. HH:MM:SS.s iff True. + """ + useSexagesimal = bool(useSexagesimal) + if useSexagesimal != self._useSexagesimal: + self._useSexagesimal = useSexagesimal + self.draw() + + def remove(self): + """Remove the AST axes and restore matplotlib's axis furniture.""" + if self._redrawTimer is not None: + self._redrawTimer.stop() + self._redrawTimer = None + for callbackId in self._callbackIds: + self._axes.callbacks.disconnect(callbackId) + self._callbackIds = [] + + self._removeArtists() + + xVisible, yVisible, spines = self._savedVisibility + self._axes.xaxis.set_visible(xVisible) + self._axes.yaxis.set_visible(yVisible) + for name, visible in spines.items(): + self._axes.spines[name].set_visible(visible) diff --git a/tests/test_display_matplotlib.py b/tests/test_display_matplotlib.py index 0fbfbc6..235f7b3 100644 --- a/tests/test_display_matplotlib.py +++ b/tests/test_display_matplotlib.py @@ -23,8 +23,334 @@ """ import unittest -import lsst.utils.tests -import lsst.afw.display as afwDisplay +import matplotlib +matplotlib.use("Agg") +import matplotlib.figure # noqa: E402 +import matplotlib.text # noqa: E402 +from matplotlib.backends.backend_agg import FigureCanvasAgg # noqa: E402 + +import lsst.utils.tests # noqa: E402 +import lsst.geom # noqa: E402 +import lsst.afw.display as afwDisplay # noqa: E402 +import lsst.afw.geom as afwGeom # noqa: E402 +import lsst.afw.image as afwImage # noqa: E402 + + +def makeTestWcs(): + """Return a rotated TAN SkyWcs for testing.""" + return afwGeom.makeSkyWcs( + crpix=lsst.geom.Point2D(100, 75), + crval=lsst.geom.SpherePoint(30.0, -45.0, lsst.geom.degrees), + cdMatrix=afwGeom.makeCdMatrix(scale=3.0*lsst.geom.arcseconds, + orientation=30*lsst.geom.degrees), + ) + + +class AstFrameSetTestCase(lsst.utils.tests.TestCase): + """Tests for the SkyWcs to pyast FrameSet conversion.""" + + def testRoundTrip(self): + """The pyast FrameSet must reproduce pixelToSky exactly.""" + from lsst.display.matplotlib.wcsAxes import astFrameSetFromWcs + + wcs = makeTestWcs() + frameSet = astFrameSetFromWcs(wcs) + for x, y in [(0.0, 0.0), (100.0, 75.0), (199.5, 149.5), (-20.0, 300.0)]: + sky = wcs.pixelToSky(x, y) + result = frameSet.tran([[x], [y]]) + self.assertEqual(result[0][0], sky[0].asRadians()) + self.assertEqual(result[1][0], sky[1].asRadians()) + + +def makeWcsAxes(**kwargs): + """Create an Agg figure with pixel-coordinate axes and a manager. + + Returns the (axes, manager) pair. + """ + from lsst.display.matplotlib.wcsAxes import WcsAxesManager, astFrameSetFromWcs + + fig = matplotlib.figure.Figure(figsize=(8, 6)) + FigureCanvasAgg(fig) + ax = fig.add_subplot(111) + ax.set_xlim(-0.5, 199.5) + ax.set_ylim(-0.5, 149.5) + wcs = makeTestWcs() + manager = WcsAxesManager(ax, astFrameSetFromWcs(wcs), **kwargs) + return ax, manager + + +class WcsAxesManagerTestCase(lsst.utils.tests.TestCase): + """Tests for WcsAxesManager drawing and artist lifecycle.""" + + def testDrawCreatesArtists(self): + ax, manager = makeWcsAxes() + self.assertGreater(len(manager.artists), 0) + # All tracked artists really are in the axes. + axesArtists = set(ax.lines) | set(ax.texts) + for artist in manager.artists: + self.assertIn(artist, axesArtists) + # Sky axis labels come from the AST SkyFrame. + texts = [a.get_text() for a in manager.artists + if isinstance(a, matplotlib.text.Text)] + self.assertIn("Right ascension", texts) + self.assertIn("Declination", texts) + + def testFurnitureHidden(self): + ax, manager = makeWcsAxes() + self.assertFalse(ax.xaxis.get_visible()) + self.assertFalse(ax.yaxis.get_visible()) + for spine in ax.spines.values(): + self.assertFalse(spine.get_visible()) + + def testTextNotClipped(self): + ax, manager = makeWcsAxes() + for artist in manager.artists: + if isinstance(artist, matplotlib.text.Text): + self.assertFalse(artist.get_clip_on()) + + def testDecimalLabelsByDefault(self): + ax, manager = makeWcsAxes() + tickLabels = [a.get_text() for a in manager.artists + if isinstance(a, matplotlib.text.Text) and + a.get_text() not in ("Right ascension", "Declination")] + self.assertGreater(len(tickLabels), 0) + for label in tickLabels: + self.assertNotIn(":", label) + # AST chooses the precision; the labels must be decimal degrees + # and adjacent labels must be distinct. + self.assertTrue(any("." in label for label in tickLabels)) + self.assertEqual(len(tickLabels), len(set(tickLabels))) + + def testSexagesimalLabels(self): + ax, manager = makeWcsAxes(useSexagesimal=True) + tickLabels = [a.get_text() for a in manager.artists + if isinstance(a, matplotlib.text.Text)] + self.assertTrue(any(":" in label for label in tickLabels)) + + def testRedrawOnLimitChange(self): + ax, manager = makeWcsAxes() + oldArtists = list(manager.artists) + ax.set_xlim(50, 150) + ax.set_ylim(40, 110) + self.assertGreater(len(manager.artists), 0) + self.assertNotEqual(manager.artists, oldArtists) + # The old artists must be gone from the axes. + axesArtists = set(ax.lines) | set(ax.texts) + for artist in oldArtists: + self.assertNotIn(artist, axesArtists) + + def testDebouncedRedraw(self): + """With a working timer, drag events coalesce into one redraw.""" + from matplotlib.backend_bases import TimerBase + + starts = [] + + class FakeTimer(TimerBase): + def _timer_start(self): + starts.append(self) + + def _timer_stop(self): + pass + + ax, manager = makeWcsAxes() + canvas = ax.get_figure().canvas + canvas.new_timer = lambda interval=None: FakeTimer(interval=interval) + + oldArtists = list(manager.artists) + for i in range(5): + ax.set_xlim(-0.5 + i, 199.5 + i) + ax.set_ylim(-0.5 + i, 149.5 + i) + # No redraw yet: the debounce timer is pending. + self.assertEqual(manager.artists, oldArtists) + self.assertGreater(len(starts), 0) + # Fire the pending timer: exactly one rebuild for the new view. + manager._redrawTimer._on_timer() + self.assertNotEqual(manager.artists, oldArtists) + axesArtists = set(ax.lines) | set(ax.texts) + for artist in oldArtists: + self.assertNotIn(artist, axesArtists) + + def testNoRedrawAfterRemove(self): + ax, manager = makeWcsAxes() + manager.remove() + ax.set_xlim(50, 150) + self.assertEqual(manager.artists, []) + self.assertEqual(len(ax.lines), 0) + + def testSetSexagesimal(self): + ax, manager = makeWcsAxes() + + def tickLabels(): + return [a.get_text() for a in manager.artists + if isinstance(a, matplotlib.text.Text)] + + self.assertFalse(any(":" in label for label in tickLabels())) + manager.setSexagesimal(True) + self.assertTrue(any(":" in label for label in tickLabels())) + manager.setSexagesimal(False) + self.assertFalse(any(":" in label for label in tickLabels())) + + def testConstructorFailureRestoresFurniture(self): + from lsst.display.matplotlib.wcsAxes import WcsAxesManager, astFrameSetFromWcs + + fig = matplotlib.figure.Figure(figsize=(8, 6)) + FigureCanvasAgg(fig) + ax = fig.add_subplot(111) + ax.set_xlim(-0.5, 199.5) + ax.set_ylim(-0.5, 149.5) + with self.assertRaises(Exception): + WcsAxesManager(ax, astFrameSetFromWcs(makeTestWcs()), + extraOptions="NoSuchAttr=1") + self.assertTrue(ax.xaxis.get_visible()) + self.assertTrue(ax.yaxis.get_visible()) + for spine in ax.spines.values(): + self.assertTrue(spine.get_visible()) + self.assertEqual(len(ax.lines) + len(ax.texts), 0) + + def testRemove(self): + ax, manager = makeWcsAxes() + manager.remove() + self.assertEqual(manager.artists, []) + self.assertEqual(len(ax.lines), 0) + self.assertEqual(len(ax.texts), 0) + self.assertTrue(ax.xaxis.get_visible()) + self.assertTrue(ax.yaxis.get_visible()) + for spine in ax.spines.values(): + self.assertTrue(spine.get_visible()) + + +class WcsAxesDisplayTestCase(lsst.utils.tests.TestCase): + """Tests for WCS axes at the afwDisplay level.""" + + frame = 100 # keep test figures separate from any other test + + def setUp(self): + afwDisplay.setDefaultBackend("matplotlib") + # A distinct frame per test so each test gets a fresh figure. + WcsAxesDisplayTestCase.frame += 1 + self.frame = WcsAxesDisplayTestCase.frame + + def tearDown(self): + afwDisplay.delAllDisplays() + + @staticmethod + def makeExposure(): + exposure = afwImage.ExposureF(200, 150) + exposure.image.array[:] = 1.0 + exposure.setWcs(makeTestWcs()) + return exposure + + def testWcsAxesByDefault(self): + display = afwDisplay.Display(frame=self.frame) + display.mtv(self.makeExposure()) + impl = display._impl + ax = impl._figure.gca() + self.assertIn(ax, impl._wcsAxesManagers) + self.assertGreater(len(impl._wcsAxesManagers[ax].artists), 0) + self.assertFalse(ax.xaxis.get_visible()) + + def testUseWcsAxesFalse(self): + display = afwDisplay.Display(frame=self.frame, useWcsAxes=False) + display.mtv(self.makeExposure()) + impl = display._impl + ax = impl._figure.gca() + self.assertEqual(impl._wcsAxesManagers, {}) + self.assertTrue(ax.xaxis.get_visible()) + + def testNoWcsMeansPixelAxes(self): + display = afwDisplay.Display(frame=self.frame) + display.mtv(afwImage.ImageF(200, 150)) + impl = display._impl + ax = impl._figure.gca() + self.assertEqual(impl._wcsAxesManagers, {}) + self.assertTrue(ax.xaxis.get_visible()) + + def testWcsGridOption(self): + display = afwDisplay.Display(frame=self.frame, wcsGrid=True) + display.mtv(self.makeExposure()) + impl = display._impl + ax = impl._figure.gca() + manager = impl._wcsAxesManagers[ax] + self.assertIn("Grid=1", manager._plotOptions()) + + def testZoomRedraws(self): + display = afwDisplay.Display(frame=self.frame) + display.mtv(self.makeExposure()) + impl = display._impl + ax = impl._figure.gca() + manager = impl._wcsAxesManagers[ax] + oldArtists = list(manager.artists) + display.zoom(4) + display.pan(100, 75) + self.assertGreater(len(manager.artists), 0) + self.assertNotEqual(manager.artists, oldArtists) + + def testErasePreservesWcsAxes(self): + display = afwDisplay.Display(frame=self.frame) + display.mtv(self.makeExposure()) + impl = display._impl + ax = impl._figure.gca() + manager = impl._wcsAxesManagers[ax] + nAstLines = sum(1 for a in manager.artists if a in set(ax.lines)) + display.dot("+", 100, 75) + self.assertGreater(len(ax.lines), nAstLines) + display.erase() + self.assertEqual(len(ax.lines), nAstLines) + axesArtists = set(ax.lines) | set(ax.texts) + for artist in manager.artists: + self.assertIn(artist, axesArtists) + + def testFallbackRestoresPixelAxes(self): + display = afwDisplay.Display(frame=self.frame, astPlotOptions="NoSuchAttr=1") + display.mtv(self.makeExposure()) + impl = display._impl + ax = impl._figure.gca() + self.assertEqual(impl._wcsAxesManagers, {}) + self.assertTrue(ax.xaxis.get_visible()) + for spine in ax.spines.values(): + self.assertTrue(spine.get_visible()) + + def testEraseRemovesPatches(self): + display = afwDisplay.Display(frame=self.frame) + display.mtv(self.makeExposure()) + ax = display._impl._figure.gca() + display.dot("o", 100, 75, size=10) + self.assertEqual(len(ax.patches), 1) + display.erase() + self.assertEqual(len(ax.patches), 0) + + def testRepeatedMtvDoesNotLeak(self): + display = afwDisplay.Display(frame=self.frame) + display.mtv(self.makeExposure()) + impl = display._impl + ax = impl._figure.gca() + nTexts = len(ax.texts) + display.mtv(self.makeExposure()) + ax = impl._figure.gca() + self.assertEqual(len(impl._wcsAxesManagers), 1) + self.assertEqual(len(ax.texts), nTexts) + + def testUseSexagesimalRedraws(self): + display = afwDisplay.Display(frame=self.frame) + display.mtv(self.makeExposure()) + impl = display._impl + ax = impl._figure.gca() + manager = impl._wcsAxesManagers[ax] + + def hasColon(): + return any(":" in a.get_text() for a in manager.artists + if isinstance(a, matplotlib.text.Text)) + + self.assertFalse(hasColon()) + display.useSexagesimal(True) + self.assertTrue(hasColon()) + + def testCloseReleasesManagers(self): + display = afwDisplay.Display(frame=self.frame) + display.mtv(self.makeExposure()) + impl = display._impl + impl._close() + self.assertEqual(impl._wcsAxesManagers, {}) class DisplayMatplotlibTestCase(unittest.TestCase):