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