Coverage for python / lsst / cell_coadds / _stitched_coadd.py: 25%
137 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-21 10:41 +0000
1# This file is part of cell_coadds.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22from __future__ import annotations
24__all__ = ("StitchedCoadd",)
26import dataclasses
27from collections.abc import Iterator, Set
28from functools import partial
29from typing import TYPE_CHECKING
31import numpy as np
33from lsst.afw.image import CoaddInputs, ExposureF, FilterLabel, PhotoCalib
34from lsst.afw.table import ExposureTable, Schema
35from lsst.geom import Box2I, Point2I
37from ._common_components import CoaddUnits, CommonComponents, CommonComponentsProperties
38from ._grid_container import GridContainer
39from ._identifiers import ObservationIdentifiers
40from ._image_planes import ImagePlanes, ViewImagePlanes
41from ._stitched_aperture_correction import StitchedApertureCorrection
42from ._stitched_exposure_catalog import StitchedExposureCatalog
43from ._stitched_image_planes import StitchedImagePlanes
44from ._stitched_psf import StitchedPsf
45from ._uniform_grid import UniformGrid
46from .typing_helpers import StitchedCoaddApCorrMap
48if TYPE_CHECKING: 48 ↛ 49line 48 didn't jump to line 49 because the condition on line 48 was never true
49 from ._multiple_cell_coadd import MultipleCellCoadd
52class StitchedCoadd(StitchedImagePlanes, CommonComponentsProperties):
53 """A lazy-evaluation coadd that stitches together images from adjacent
54 cells.
56 Parameters
57 ----------
58 cell_coadd : `MultipleCellCoadd`
59 Cell-based coadd to stitch together.
60 bbox : `Box2I`, optional
61 The region over which a contiguous coadd is desired. Defaults to
62 ``cell_coadd.outer_bbox``.
64 Notes
65 -----
66 This class simply inserts subimages from each cell into the full image,
67 doing so when an attribute is first accessed to avoid stitching together
68 planes that may never be accessed.
70 A `StitchedCoadd` cannot be serialized in FITS format directly. Instead,
71 the recommended way is to serialize the `MultipleCellCoadd` instance that
72 was used to construct the object and reconstruct the `StitchedCoadd` by
73 calling the `stitch` method on it. A less recommended way is to call the
74 `asExposure` method to get an `lsst.afw.image.Exposure` object and persist
75 that to the disk.
76 """
78 def __init__(self, cell_coadd: MultipleCellCoadd, *, bbox: Box2I | None = None):
79 super().__init__()
80 if bbox is None:
81 bbox = cell_coadd.outer_bbox
82 elif not cell_coadd.outer_bbox.contains(bbox):
83 raise ValueError(
84 f"Cell coadd inner bounding box {cell_coadd.outer_bbox} does not "
85 f"contain stitch target area {bbox}."
86 )
87 self._bbox = bbox
88 self._cell_coadd = cell_coadd
89 self._psf: StitchedPsf | None = None
90 self._ap_corr_map: StitchedCoaddApCorrMap | None = None
91 self._ccds: StitchedExposureCatalog | None = None
92 self._visits: StitchedExposureCatalog | None = None
93 self._common = cell_coadd.common
95 @property
96 def bbox(self) -> Box2I:
97 # Docstring inherited.
98 return self._bbox
100 @property
101 def grid(self) -> UniformGrid:
102 """Object that defines the piecewise grid (of inner cell regions) that
103 this object stitches together.
105 This may include cells outside the region covered by these image
106 planes.
107 """
108 return self._cell_coadd.grid
110 def _iter_cell_planes(self) -> Iterator[ImagePlanes]:
111 # Docstring inherited.
112 x_max, y_max = self.grid.index(self.bbox.getMax())
114 for cell in self._cell_coadd.cells.values():
115 bbox = cell.inner.bbox
116 if cell.identifiers.cell.x == 0:
117 # This is a special case for the first column of cells.
118 bbox.include(Point2I(cell.outer.bbox.beginX, cell.outer.bbox.centerY))
119 elif cell.identifiers.cell.x == x_max:
120 # This is a special case for the last column of cells.
121 bbox.include(Point2I(cell.outer.bbox.endX, cell.outer.bbox.centerY))
123 if cell.identifiers.cell.y == 0:
124 # This is a special case for the last row of cells.
125 bbox.include(Point2I(cell.outer.bbox.centerX, cell.outer.bbox.beginY))
126 elif cell.identifiers.cell.y == y_max:
127 # This is a special case for the first row of cells.
128 bbox.include(Point2I(cell.outer.bbox.centerX, cell.outer.bbox.endY))
130 bbox.clip(cell.outer.bbox)
131 make_view = partial(cell.make_view, bbox=bbox)
132 yield ViewImagePlanes(cell.outer, bbox=bbox, make_view=make_view)
134 @property
135 def n_noise_realizations(self) -> int:
136 # Docstring inherited.
137 return self._cell_coadd.n_noise_realizations
139 @property
140 def mask_fraction_names(self) -> Set[str]:
141 # Docstring inherited.
142 return self._cell_coadd.mask_fraction_names
144 @property
145 def psf(self) -> StitchedPsf:
146 """The piecewise PSF of this image."""
147 if self._psf is None:
148 self._psf = StitchedPsf(
149 self._cell_coadd.cells.rebuild_transformed(lambda cell: cell.psf_image.convertD()),
150 self._cell_coadd.grid,
151 )
152 return self._psf
154 @property
155 def common(self) -> CommonComponents:
156 # Docstring inherited.
157 return self._cell_coadd.common
159 def set_cell_edges(self, *, edge_width: int = 1, edge_mask_name: str = "CELL_EDGE") -> None:
160 """Set a mask bit indicating the inner cell edges.
162 Parameters
163 ----------
164 edge_width : `int`, optional
165 The width of the edge region to flag, in pixels.
166 edge_mask_name : `str`, optional
167 The name of the mask plane to add for the edge region.
168 """
169 for cell in self._cell_coadd.cells.values():
170 cell._set_cell_edges(edge_width=edge_width, edge_mask_name=edge_mask_name)
172 def _make_coadd_inputs(self) -> CoaddInputs:
173 """Make a coadd inputs object for this stitched coadd.
175 Returns
176 -------
177 coadd_inputs : `lsst.afw.image.CoaddInputs`
178 The coadd inputs for this stitched coadd.
179 """
180 visit_schema = ExposureTable.makeMinimalSchema()
181 visit_schema.addField("visit", type=np.int64, doc="Visit ID")
183 ccd_schema = Schema(visit_schema)
184 ccd_schema.addField("ccd", type=np.int32, doc="Detector ID")
185 ccd_schema.addField("weight", type=np.float32, doc="Weight for this (visit, detector) in the coadd.")
187 coadd_inputs = CoaddInputs(visit_schema, ccd_schema)
189 for obs_id in self.visits:
190 visit_record = coadd_inputs.visits.addNew()
191 visit_record["visit"] = obs_id.visit
192 visit_record.setBBox(Box2I()) # Dummy bbox; not used.
194 # Pre-process the weights to avoid redundant lookups and matching.
195 weights: dict[ObservationIdentifiers, float] = {}
196 for cell in self._cell_coadd.cells.values():
197 for obs_id, coadd_input in cell.inputs.items():
198 weights[obs_id] = coadd_input.weight
199 if len(weights) == len(self._cell_coadd.common.visit_polygons):
200 break
202 for obs_id, polygon in self._cell_coadd.common.visit_polygons.items():
203 ccd_record = coadd_inputs.ccds.addNew()
204 ccd_record["visit"] = obs_id.visit
205 ccd_record["ccd"] = obs_id.detector
206 # ObservationIdentifiers in the inputs is a proper subset of those
207 # denoting the visit polygons. For edgeless coadds, some
208 # identifiers may not be present in weights and we take that to
209 # mean that they were excluded from the coadds (and hence have
210 # zero weight).
211 ccd_record["weight"] = weights.get(obs_id, 0.0)
212 ccd_record.setBBox(Box2I()) # Dummy bbox; not used.
213 ccd_record.validPolygon = polygon
215 return coadd_inputs
217 def asExposure(
218 self,
219 *,
220 noise_index: int | None = None,
221 ) -> ExposureF:
222 """Return an `lsst.afw.image.Exposure` view of this piecewise image.
224 Parameters
225 ----------
226 noise_index : `int` or `None`, optional
227 If specified, return an exposure containing the specified noise
228 realization instead of the data image.
230 Returns
231 -------
232 exposure : `lsst.afw.image.ExposureF`
233 The stitched exposure.
235 Notes
236 -----
237 The returned exposure will have its coaddInputs set. However, its
238 internal representation is very different from those that are genuine
239 afw objects. The records in the `ccds` and `visits` tables do not
240 contain the individual bounding boxes and WCSs. Instead, the
241 validPolygon is already in the coadd coordinates. To query the subset
242 of inputs, use the `subset_containing_ccds` and
243 `subset_containing_visits` methods on the returned exposure's
244 `getInfo().getCoaddInputs()` object.
245 """
246 result = ExposureF(self.asMaskedImage(noise_index=noise_index))
247 result.getInfo().setCoaddInputs(self._make_coadd_inputs())
248 # Exposure components derived from "common" components are all simple.
249 result.setWcs(self._cell_coadd.wcs)
250 result.setFilter(FilterLabel(band=self.band))
251 if self.units is CoaddUnits.nJy:
252 result.setPhotoCalib(PhotoCalib(1.0))
253 result.metadata["BUNIT"] = "nJy"
255 # We can't do result.setId here, because:
256 #
257 # - we don't know whether this should be a packed tract+patch+band ID
258 # or just a tract+patch ID;
259 #
260 # - we don't know how to pack the information we have anyway.
261 #
262 # Maybe DM-31924 will provide a solution to at least the latter.
263 # result.setId(self._cell_coadd.identifiers.patch)
265 # We could add CoaddInputs here, but without WCS, PSF, etc in them;
266 # it's not clear that's good enough or even useful, given that the cell
267 # provide a much more clean view of what the inputs are at any given
268 # point.
270 # PSF is the first of many components that need piecewise
271 # implementations. More to do here for at least aperture corrections
272 # and transmission curves.
273 result.setPsf(self.psf)
275 return result
277 @property
278 def ap_corr_map(self) -> StitchedCoaddApCorrMap:
279 """Stitch the aperture correction maps from the cell coadd.
281 This converts the aperture correction maps from each cell into a single
282 `ApCorrMap` that quacks like `lsst.afw.image.ApCorrMap`. The resulting
283 object has the fields to correct as keys and a
284 `StitchedApertureCorrection` for each field.
286 Notes
287 -----
288 These cannot be attached to an `~lsst.afw.image.Exposure` object.
289 """
290 if self._ap_corr_map is None:
291 ap_corr_map: dict[str, StitchedApertureCorrection] = {}
292 field_names = self._cell_coadd.cells.arbitrary.aperture_corrected_algorithms
293 for field_name in field_names:
294 gc = GridContainer[float](shape=self.grid.shape)
295 for scc in self._cell_coadd.cells.values():
296 gc[scc.identifiers.cell] = scc.aperture_correction_map[field_name]
297 ap_corr_map[field_name] = StitchedApertureCorrection(self.grid, gc)
299 self._ap_corr_map = ap_corr_map
301 return self._ap_corr_map
303 @property
304 def ccds(self) -> StitchedExposureCatalog:
305 """An ExposureCatalog-like object that contains the information about
306 the ``ObservationIdentifiers`` at any point on the coadd. It can also
307 be used as an iterable to iterate over the unique (visit, detector)
308 combinations.
309 """
310 if self._ccds is None:
311 gc = GridContainer[tuple[ObservationIdentifiers, ...]](shape=self.grid.shape)
312 for idx, scc in self._cell_coadd.cells.items():
313 gc[idx] = tuple(scc.inputs.keys())
314 self._ccds = StitchedExposureCatalog(self.grid, gc)
316 return self._ccds
318 @property
319 def visits(self) -> StitchedExposureCatalog:
320 """An ExposureCatalog-like object that contains the information about
321 the ``ObservationIdentifiers`` at any point on the coadd. It can also
322 be used as an iterable to iterate over the unique visits.
323 """
324 if self._visits is None:
325 gc = GridContainer[tuple[ObservationIdentifiers, ...]](shape=self.grid.shape)
326 for idx, scc in self._cell_coadd.cells.items():
327 gc[idx] = tuple(
328 sorted({dataclasses.replace(ccd_input, detector=-1) for ccd_input in scc.inputs})
329 )
330 self._visits = StitchedExposureCatalog(self.grid, gc)
332 return self._visits