Coverage for python/lsst/cell_coadds/_stitched_coadd.py: 25%
137 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-20 02:20 -0800
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 ccd_record["weight"] = weights[obs_id]
207 ccd_record.setBBox(Box2I()) # Dummy bbox; not used.
208 ccd_record.validPolygon = polygon
210 return coadd_inputs
212 def asExposure(
213 self,
214 *,
215 noise_index: int | None = None,
216 ) -> ExposureF:
217 """Return an `lsst.afw.image.Exposure` view of this piecewise image.
219 Parameters
220 ----------
221 noise_index : `int` or `None`, optional
222 If specified, return an exposure containing the specified noise
223 realization instead of the data image.
225 Returns
226 -------
227 exposure : `lsst.afw.image.ExposureF`
228 The stitched exposure.
230 Notes
231 -----
232 The returned exposure will have its coaddInputs set. However, its
233 internal representation is very different from those that are genuine
234 afw objects. The records in the `ccds` and `visits` tables do not
235 contain the individual bounding boxes and WCSs. Instead, the
236 validPolygon is already in the coadd coordinates. To query the subset
237 of inputs, use the `subset_containing_ccds` and
238 `subset_containing_visits` methods on the returned exposure's
239 `getInfo().getCoaddInputs()` object.
240 """
241 result = ExposureF(self.asMaskedImage(noise_index=noise_index))
242 result.getInfo().setCoaddInputs(self._make_coadd_inputs())
243 # Exposure components derived from "common" components are all simple.
244 result.setWcs(self._cell_coadd.wcs)
245 result.setFilter(FilterLabel(band=self.band))
246 if self.units is CoaddUnits.nJy:
247 result.setPhotoCalib(PhotoCalib(1.0))
248 result.metadata["BUNIT"] = "nJy"
250 # We can't do result.setId here, because:
251 #
252 # - we don't know whether this should be a packed tract+patch+band ID
253 # or just a tract+patch ID;
254 #
255 # - we don't know how to pack the information we have anyway.
256 #
257 # Maybe DM-31924 will provide a solution to at least the latter.
258 # result.setId(self._cell_coadd.identifiers.patch)
260 # We could add CoaddInputs here, but without WCS, PSF, etc in them;
261 # it's not clear that's good enough or even useful, given that the cell
262 # provide a much more clean view of what the inputs are at any given
263 # point.
265 # PSF is the first of many components that need piecewise
266 # implementations. More to do here for at least aperture corrections
267 # and transmission curves.
268 result.setPsf(self.psf)
270 return result
272 @property
273 def ap_corr_map(self) -> StitchedCoaddApCorrMap:
274 """Stitch the aperture correction maps from the cell coadd.
276 This converts the aperture correction maps from each cell into a single
277 `ApCorrMap` that quacks like `lsst.afw.image.ApCorrMap`. The resulting
278 object has the fields to correct as keys and a
279 `StitchedApertureCorrection` for each field.
281 Notes
282 -----
283 These cannot be attached to an `~lsst.afw.image.Exposure` object.
284 """
285 if self._ap_corr_map is None:
286 ap_corr_map: dict[str, StitchedApertureCorrection] = {}
287 field_names = self._cell_coadd.cells.arbitrary.aperture_corrected_algorithms
288 for field_name in field_names:
289 gc = GridContainer[float](shape=self.grid.shape)
290 for scc in self._cell_coadd.cells.values():
291 gc[scc.identifiers.cell] = scc.aperture_correction_map[field_name]
292 ap_corr_map[field_name] = StitchedApertureCorrection(self.grid, gc)
294 self._ap_corr_map = ap_corr_map
296 return self._ap_corr_map
298 @property
299 def ccds(self) -> StitchedExposureCatalog:
300 """An ExposureCatalog-like object that contains the information about
301 the ``ObservationIdentifiers`` at any point on the coadd. It can also
302 be used as an iterable to iterate over the unique (visit, detector)
303 combinations.
304 """
305 if self._ccds is None:
306 gc = GridContainer[tuple[ObservationIdentifiers, ...]](shape=self.grid.shape)
307 for idx, scc in self._cell_coadd.cells.items():
308 gc[idx] = tuple(scc.inputs.keys())
309 self._ccds = StitchedExposureCatalog(self.grid, gc)
311 return self._ccds
313 @property
314 def visits(self) -> StitchedExposureCatalog:
315 """An ExposureCatalog-like object that contains the information about
316 the ``ObservationIdentifiers`` at any point on the coadd. It can also
317 be used as an iterable to iterate over the unique visits.
318 """
319 if self._visits is None:
320 gc = GridContainer[tuple[ObservationIdentifiers, ...]](shape=self.grid.shape)
321 for idx, scc in self._cell_coadd.cells.items():
322 gc[idx] = tuple(
323 sorted({dataclasses.replace(ccd_input, detector=-1) for ccd_input in scc.inputs})
324 )
325 self._visits = StitchedExposureCatalog(self.grid, gc)
327 return self._visits