Coverage for python / lsst / analysis / tools / tasks / makeTractPostageStamp.py: 21%
88 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-05 18:53 +0000
1# This file is part of analysis_tools.
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__all__ = (
23 "MakeTractPostageStampConfig",
24 "MakeTractPostageStampTask",
25)
27import copy
29import lsst.pipe.base as pipeBase
30import matplotlib.cm as cm
31import matplotlib.patheffects as pathEffects
32import numpy as np
33from lsst.pipe.base import connectionTypes as ct
34from lsst.skymap import BaseSkyMap
35from lsst.utils.plotting import make_figure
37from ..utils import getPatchCorners, getTractCorners
40class MakeTractPostageStampConnections(
41 pipeBase.PipelineTaskConnections,
42 dimensions=("skymap", "tract", "band"),
43 defaultTemplates={
44 "inputName": "deepCoadd_calexpBin",
45 "outputName": "deepTract_PostageStamp",
46 },
47):
48 data = ct.Input(
49 doc="Binned deepCoadd calibrated exposures to read from the butler.",
50 name="{inputName}",
51 storageClass="ExposureF",
52 deferLoad=True,
53 dimensions=(
54 "skymap",
55 "tract",
56 "patch",
57 "band",
58 ),
59 multiple=True,
60 )
62 skymap = ct.Input(
63 doc="The skymap that covers the tract that the data is from.",
64 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
65 storageClass="SkyMap",
66 dimensions=("skymap",),
67 )
69 postageStamp = ct.Output(
70 doc="A postagestamp composite image of the tract.",
71 name="{outputName}",
72 storageClass="Plot",
73 dimensions=("skymap", "tract"),
74 )
77class MakeTractPostageStampConfig(
78 pipeBase.PipelineTaskConfig, pipelineConnections=MakeTractPostageStampConnections
79):
80 pass
83class MakeTractPostageStampTask(pipeBase.PipelineTask):
85 ConfigClass = MakeTractPostageStampConfig
86 _DefaultName = "makeTractPostageStamp"
88 def runQuantum(self, butlerQC, inputRefs, outputRefs):
89 """Takes a set of coadded patch Exposures and displays them
90 in their corresponding positions within a tract.
91 Empty patches - those that do not have any coverage - are shown
92 as hatched squares.
94 Parameters
95 ----------
96 butlerQC : `lsst.pipe.base.QuantumContext`
97 A butler which is specialized to operate in the context of a
98 `lsst.daf.butler.Quantum`.
99 inputRefs : `lsst.pipe.base.InputQuantizedConnection`
100 Datastructure containing named attributes 'data and 'skymap'.
101 The values of these attributes are the corresponding
102 `lsst.daf.butler.DatasetRef` objects defined in the corresponding
103 `PipelineTaskConnections` class.
104 outputRefs : `lsst.pipe.base.OutputQuantizedConnection`
105 Datastructure containing named attribute 'postageStamp'.
106 The value of this attribute is the corresponding
107 `lsst.daf.butler.DatasetRef` object defined in the corresponding
108 `PipelineTaskConnections` class.
109 """
111 inputs = butlerQC.get(inputRefs)
112 patches = inputs["data"]
113 skymap = inputs["skymap"]
114 band = butlerQC.quantum.dataId["band"]
115 tract = butlerQC.quantum.dataId["tract"]
117 fig = self.makeTractPostageStamp(skymap, tract, patches, band)
119 butlerQC.put(pipeBase.Struct(postageStamp=fig), outputRefs)
121 def makeTractPostageStamp(self, skymap, tract, patches, band):
122 """Takes the coadd patch exposures and displays them on a
123 set of axes. The axes boundaries are those of the tract.
124 Patches are annoted with their patch identification number.
125 Empty patches - those that do not have any coverage - are shown
126 as hatched squares.
128 Parameters
129 ----------
130 skymap : `lsst.skymap`
131 tract : `int`
132 Identification number of tract.
133 patches : `list` [`DeferredDatasetHandle`]
134 List of handles for patch coadd exposures to display.
135 band : `str`
136 Filter band. Only used to annotate the plot.
138 Returns
139 -------
140 fig : `matplotlib.figure.Figure`
141 Plot displaying the tract with coadd exposures displayed.
142 """
144 tractInfo = skymap.generateTract(tract)
145 tractCorners = getTractCorners(skymap, tract)
146 tractRas = [ra for (ra, dec) in tractCorners]
147 RaSpansZero = max(tractRas) > 360.0
149 cmap = cm.grey
150 cmap.set_bad(alpha=0)
151 cmapred = cm.Reds
152 cmapred.set_bad(alpha=0)
154 fig = make_figure(dpi=300)
155 ax = fig.add_subplot(111)
157 emptyPatches = np.arange(tractInfo.getNumPatches()[0] * tractInfo.getNumPatches()[1]).tolist()
159 for patch in patches:
161 patchId = patch.dataId["patch"]
162 emptyPatches.remove(patchId)
164 # Create the patch axes at the appropriate location in tract:
165 patchCorners = getPatchCorners(tractInfo, patchId)
166 ras = [ra for (ra, dec) in patchCorners]
167 decs = [dec for (ra, dec) in patchCorners]
169 # Account for the RA wrapping using negative RA values.
170 # This is rectified when the final axes are built.
171 if RaSpansZero:
172 ras = [ra - 360 if ra > 180.0 else ra for ra in ras]
173 Extent = (max(ras), min(ras), max(decs), min(decs))
174 ax.plot(
175 [min(ras), max(ras), max(ras), min(ras), min(ras)],
176 [min(decs), min(decs), max(decs), max(decs), min(decs)],
177 "k",
178 lw=0.5,
179 )
181 # Fetch the images and display:
182 im = patch.get(component="image").array
183 mask = patch.get(component="mask").array & 2**8 > 0
185 imdata = copy.deepcopy(im)
186 imdata[mask] = np.nan
187 med = np.nanmedian(imdata)
188 mad = np.nanmedian(np.fabs(im - med))
189 vmin = med - 1 * mad
190 vmax = med + 15 * mad
191 ax.imshow(imdata, extent=Extent, vmin=vmin, vmax=vmax, cmap=cmap)
193 # Highlight regions of "NODATA".
194 imnodata = copy.deepcopy(im)
195 imnodata[~mask] = np.nan
196 ax.imshow(imnodata, extent=Extent, alpha=0.8, cmap=cmapred)
198 ax.annotate(
199 patchId,
200 (np.mean(ras), np.mean(decs)),
201 color="k",
202 ha="center",
203 va="center",
204 path_effects=[pathEffects.withStroke(linewidth=2, foreground="w")],
205 )
207 # Indicate blank patches as hatched regions
208 for patchId in emptyPatches:
210 patchCorners = getPatchCorners(tractInfo, patchId)
211 ras = [ra for (ra, dec) in patchCorners]
212 decs = [dec for (ra, dec) in patchCorners]
214 # Account for the RA wrapping using negative RA values.
215 if RaSpansZero:
216 ras = [ra - 360 if ra > 180.0 else ra for ra in ras]
218 Extent = (max(ras), min(ras), max(decs), min(decs))
219 ax.plot(
220 [min(ras), max(ras), max(ras), min(ras), min(ras)],
221 [min(decs), min(decs), max(decs), max(decs), min(decs)],
222 "k",
223 lw=0.5,
224 )
226 cs = ax.contourf(np.ones((10, 10)), 1, hatches=["xx"], extent=Extent, colors="none")
227 cs.set_edgecolors("red")
228 ax.annotate(
229 patchId,
230 (np.mean(ras), np.mean(decs)),
231 color="k",
232 ha="center",
233 va="center",
234 path_effects=[pathEffects.withStroke(linewidth=2, foreground="w")],
235 )
237 # Draw axes around the entire tract:
238 ax.set_xlabel("R.A. (deg)", fontsize=20)
239 ax.set_ylabel("Dec. (deg)", fontsize=20)
241 tractRas = [ra for (ra, dec) in tractCorners]
242 # Account for the RA wrapping using negative RA values.
243 if RaSpansZero:
244 tractRas = [ra - 360.0 for ra in tractRas]
246 ax.set_xlim(max(tractRas), min(tractRas))
247 ticks = [t for t in ax.get_xticks() if t >= min(tractRas) and t <= max(tractRas)]
249 # Rectify potential negative RA values via tick labels
250 tickLabels = [f"{t % 360:.1f}" for t in ticks]
251 ax.set_xticks(ticks, tickLabels)
253 tractDecs = [dec for (ra, dec) in tractCorners]
254 ax.set_ylim(min(tractDecs), max(tractDecs))
256 ax.tick_params(axis="both", labelsize=10, length=0, pad=1.5)
257 ax.set_title(str(tract) + ": " + band, fontsize=20)
258 fig.canvas.draw()
260 return fig