Coverage for python/lsst/images/cells/_aperture_corrections.py: 54%
197 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-11-08 01:42 -0800
1# This file is part of lsst-images.
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# Use of this source code is governed by a 3-clause BSD-style
10# license that can be found in the LICENSE file.
12from __future__ import annotations
14__all__ = (
15 "CellApertureCorrectionMapSerializationModel",
16 "CellField",
17 "approximate_aperture_corrections_from_legacy",
18)
20from collections.abc import Mapping
21from typing import TYPE_CHECKING, Any, ClassVar, final
23import astropy.table
24import astropy.units
25import numpy as np
26import pydantic
28from .._cell_grid import CellGridBounds, CellIJ
29from .._geom import BoundsError, Box
30from .._image import Image
31from ..fields import BaseField
32from ..serialization import (
33 ArchiveReadError,
34 ArchiveTree,
35 InputArchive,
36 InvalidParameterError,
37 OutputArchive,
38 TableModel,
39)
41if TYPE_CHECKING:
42 try:
43 from lsst.afw.math import BoundedField as LegacyBoundedField
44 from lsst.cell_coadds import StitchedApertureCorrection as LegacyStitchedApertureCorrection
45 except ImportError:
46 type LegacyBoundedField = Any # type: ignore[no-redef]
47 type LegacyStitchedApertureCorrection = Any # type: ignore[no-redef]
50@final
51class CellField(BaseField):
52 """A piecewise 2-d function on a cell-coadd grid.
54 Parameters
55 ----------
56 bounds
57 Description of the cell grid and any missing cells. Array entries for
58 missing cells should be NaN.
59 array
60 A 2-d array of cell values with shape
61 ``bounds.subgrid_size.as_tuple()``.
62 unit
63 Units of the field values, or `None` if dimensionless.
65 Notes
66 -----
67 `CellField` is not directly serializable and is not included in the
68 ``Field`` union type alias as a result. A `~collections.abc.Mapping` of
69 `CellField` is instead serializable via
70 `CellApertureCorrectionMapSerializationModel`.
71 """
73 def __init__(
74 self, bounds: CellGridBounds, array: np.ndarray, unit: astropy.units.UnitBase | None = None
75 ) -> None:
76 self._array = array
77 self._bounds = bounds
78 self._unit = unit
79 if self._array.shape != self._bounds.subgrid_size.as_tuple(): 79 ↛ 80line 79 didn't jump to line 80 because the condition on line 79 was never true
80 raise ValueError(
81 f"Array shape ({self._array.shape}) differs from subgrid size ({self._bounds.subgrid_size})."
82 )
84 __hash__ = None # type: ignore[assignment]
86 @property
87 def bounds(self) -> CellGridBounds:
88 return self._bounds
90 @property
91 def unit(self) -> astropy.units.UnitBase | None:
92 return self._unit
94 @property
95 def is_constant(self) -> bool:
96 indices = iter(self._bounds.cell_indices())
97 try:
98 first = self.value_in_cell(next(indices))
99 except StopIteration:
100 return True
101 for other_index in indices:
102 if self.value_in_cell(other_index) != first:
103 return False
104 return True
106 def value_in_cell(self, key: CellIJ) -> float:
107 """Return the value of the field in the cell with the given index.
109 Parameters
110 ----------
111 key
112 Index of the cell to evaluate.
113 """
114 if key in self._bounds.missing:
115 raise BoundsError(f"Cell {key} is missing for this field.")
116 if not self._bounds.contains_cell(key):
117 raise BoundsError(f"Cell {key} is out of bounds for this field.")
118 index = key - self._bounds.subgrid_start
119 try:
120 return self._array[index.i, index.j]
121 except IndexError:
122 raise BoundsError(f"Cell {key} is out of bounds for this field.") from None
124 def quantity_in_cell(self, key: CellIJ) -> astropy.units.Quantity:
125 """Return the quantity (value with units) of the field in the cell
126 with the given index.
128 Parameters
129 ----------
130 key
131 Index of the cell to evaluate.
132 """
133 return astropy.units.Quantity(self.value_in_cell(key), self._unit)
135 def _evaluate(
136 self, *, x: np.ndarray, y: np.ndarray, quantity: bool
137 ) -> np.ndarray | astropy.units.Quantity:
138 # This implementation is optimized for the case where there are many
139 # more evaluation points than cells. We could switch to an
140 # implementation that zip-broadcast-iterates over x and y when that is
141 # not the case, but that feels like a premature optimization right now.
142 result = np.full(np.broadcast_shapes(y.shape, x.shape), np.nan, dtype=np.float64)
143 for cell_index in self._bounds.cell_indices():
144 cell_bbox = self._bounds.grid.bbox_of(cell_index)
145 result[cell_bbox.contains(x=x, y=y)] = self.value_in_cell(cell_index)
146 if quantity:
147 return astropy.units.Quantity(result, self._unit)
148 return result
150 def render(self, bbox: Box | None = None, *, dtype: np.typing.DTypeLike | None = None) -> Image:
151 if bbox is None:
152 bbox = self._bounds.bbox
153 bounds = self._bounds
154 else:
155 bounds = self._bounds[bbox]
156 result = Image(np.nan, bbox=bbox, dtype=dtype, unit=self._unit)
157 for cell_index in bounds.cell_indices():
158 cell_bbox = self._bounds.grid.bbox_of(cell_index).intersection(bbox)
159 result[cell_bbox].array = self.value_in_cell(cell_index)
160 return result
162 def _multiply_constant(
163 self, factor: float | astropy.units.Quantity | astropy.units.UnitBase
164 ) -> CellField:
165 factor, unit = self._handle_factor_units(factor)
166 return CellField(self._bounds, self._array * factor, unit=unit)
168 @staticmethod
169 def from_legacy_aperture_correction(
170 legacy: LegacyStitchedApertureCorrection, bounds: CellGridBounds
171 ) -> CellField:
172 """Convert from a legacy `lsst.cell_coadds.StitchedApertureCorrection`.
174 Parameters
175 ----------
176 legacy
177 Legacy field to convert.
178 bounds
179 The grid and bounds of the returned field.
180 """
181 array = np.full(bounds.subgrid_size.as_tuple(), np.nan, dtype=np.float64)
182 new_missing: set[CellIJ] = set()
183 for cell_index in bounds.cell_indices():
184 array_index = cell_index - bounds.subgrid_start
185 try:
186 value = legacy.gc[cell_index.to_legacy()]
187 except KeyError:
188 value = np.nan
189 new_missing.add(cell_index)
190 array[array_index.i, array_index.j] = value
191 if new_missing:
192 bounds = CellGridBounds(
193 grid=bounds.grid, missing=frozenset(bounds.missing | new_missing), bbox=bounds.bbox
194 )
195 return CellField(bounds, array)
197 def to_legacy_aperture_correction(self) -> LegacyStitchedApertureCorrection:
198 """Convert to a legacy
199 `lsst.cell_coadds.StitchedApertureCorrection`.
200 """
201 from lsst.cell_coadds import GridContainer, StitchedApertureCorrection
203 grid = self.bounds.grid.to_legacy()
204 gc = GridContainer[float](grid.shape)
205 for cell_index in self.bounds.cell_indices():
206 gc[cell_index.to_legacy()] = self.value_in_cell(cell_index)
207 return StitchedApertureCorrection(grid, gc)
209 @staticmethod
210 def approximate_from_legacy(
211 legacy: LegacyBoundedField, bounds: CellGridBounds, default_is_missing: bool = True
212 ) -> CellField:
213 """Convert from a legacy bounded field by evaluating it at cell
214 centers and assuming variation within a cell is approximately
215 constant.
217 When the center sample is unusable (e.g. the center falls outside
218 the field's bounding box), the field is instead sampled at the
219 cell's four corners and the usable samples are averaged. Cells
220 with no usable sample at all are marked missing.
222 Parameters
223 ----------
224 legacy
225 Legacy bounded field to approximate.
226 bounds
227 Cell grid and known-missing cells.
228 default_is_missing
229 If `True` and ``legacy.getThrowOnMissing()`` returns `False`,
230 treat values that are identical to ``legacy.getDefault()`` as
231 missing. This covers some confusing
232 `lsst.meas.algorithms.CoaddBoundedField` behavior for backwards
233 compatibility.
234 """
235 from lsst.pex.exceptions import DomainError
237 get_throw_on_missing = getattr(legacy, "getThrowOnMissing", None)
238 get_default = getattr(legacy, "getDefault", None)
239 default = (
240 get_default()
241 if default_is_missing
242 and get_throw_on_missing is not None
243 and get_default is not None
244 and not get_throw_on_missing()
245 else None
246 )
248 def usable(value: float) -> bool:
249 return np.isfinite(value) and (default is None or value != default)
251 array = np.full(bounds.subgrid_size.as_tuple(), np.nan, dtype=np.float64)
252 missing: set[CellIJ] = set()
253 for cell_index in bounds.cell_indices():
254 cell_bbox = bounds.grid.bbox_of(cell_index)
255 try:
256 value = legacy.evaluate(cell_bbox.center.to_legacy_float_point())
257 except DomainError:
258 value = np.nan
259 if not usable(value):
260 corner_values: list[float] = []
261 for corner in cell_bbox.boundary():
262 try:
263 corner_value = legacy.evaluate(corner.to_legacy_float_point())
264 except DomainError:
265 continue
266 if not usable(corner_value):
267 continue
268 corner_values.append(corner_value)
269 if not corner_values:
270 missing.add(cell_index)
271 continue
272 else:
273 value = np.mean(corner_values)
274 array_index = cell_index - bounds.subgrid_start
275 array[array_index.i, array_index.j] = value
276 if missing:
277 bounds = CellGridBounds(
278 grid=bounds.grid, missing=frozenset(bounds.missing | missing), bbox=bounds.bbox
279 )
280 return CellField(bounds, array)
283def approximate_aperture_corrections_from_legacy(
284 legacy_map: Mapping[str, LegacyBoundedField],
285 bounds: CellGridBounds,
286 flux_suffix: str = "_instFlux",
287 err_suffix: str = "_instFluxErr",
288) -> dict[str, CellField]:
289 """Convert a mapping of legacy aperture-correction bounded fields to
290 `CellField` instances by approximating them on a cell grid.
292 Parameters
293 ----------
294 legacy_map
295 Mapping of aperture-correction names to legacy bounded fields, as
296 returned by a legacy `lsst.afw.image.Exposure`'s ``getApCorrMap()``.
297 bounds
298 Cell grid and known-missing cells to convert onto.
299 flux_suffix
300 Suffix of the flux member of an aperture-correction name pair.
301 err_suffix
302 Suffix of the error member of an aperture-correction name pair.
304 Returns
305 -------
306 `dict` [`str`, `CellField`]
307 Mapping with the same keys, holding the converted fields.
309 Notes
310 -----
311 This routine relies on the fact that aperture correction fields come in
312 ``*_instFlux`` + ``*_instFluxErr`` pairs to resolve the fact that zero is
313 a possible valid (albeit suspicious) value for the latter *and* the
314 sentinel used for missing data. We assume a zero ``*_instFlux``
315 evaluation represents bad data, but a zero ``*_instFluxErr`` is bad data
316 only if the corresponding ``*_instFlux`` evaluation is also bad data.
317 """
318 converted: dict[str, CellField] = {}
319 for name, legacy_field in legacy_map.items():
320 if name.endswith(err_suffix):
321 continue
322 converted[name] = CellField.approximate_from_legacy(legacy_field, bounds)
323 for name, legacy_field in legacy_map.items():
324 if not name.endswith(err_suffix):
325 continue
326 partner = converted.get(name[: -len(err_suffix)] + flux_suffix)
327 err_bounds = partner.bounds if partner is not None else bounds
328 converted[name] = CellField.approximate_from_legacy(
329 legacy_field, err_bounds, default_is_missing=False
330 )
331 return converted
334class CellApertureCorrectionMapSerializationModel(ArchiveTree):
335 """A serialization model for a `~collections.abc.Mapping` of `CellField`,
336 which is used to represent aperture corrections for cell-based coadds.
337 """
339 SCHEMA_NAME: ClassVar[str] = "cell_aperture_correction_map"
340 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
341 MIN_READ_VERSION: ClassVar[int] = 1
342 PUBLIC_TYPE: ClassVar[type] = dict
344 table: TableModel = pydantic.Field(
345 description="Table with one row for each cell and different photometry algorithms in columns."
346 )
347 bounds: CellGridBounds = pydantic.Field(
348 description=(
349 "Description of the cell grid and any missing cells. Array entries for "
350 "missing cells should be NaN."
351 ),
352 )
354 @staticmethod
355 def serialize(
356 aperture_correction_map: Mapping[str, CellField], archive: OutputArchive[Any]
357 ) -> CellApertureCorrectionMapSerializationModel | None:
358 if not aperture_correction_map: 358 ↛ 359line 358 didn't jump to line 359 because the condition on line 358 was never true
359 return None
360 bounds = next(iter(aperture_correction_map.values())).bounds
361 for field in aperture_correction_map.values():
362 if field.bounds.grid != bounds.grid:
363 raise ValueError("Cell aperture corrections do not have a consistent grid.")
364 if field.bounds.bbox != bounds.bbox:
365 raise ValueError("Cell aperture corrections do not have a consistent bounding box.")
366 if not all(field.bounds == bounds for field in aperture_correction_map.values()):
367 # In rare cases some field can have more missing cells from the
368 # others. This isn't worth a full denormalization into full
369 # per-field bounds (especially since that's a schema change).
370 # Instead, the shared, serialized 'bounds' only marks a cell as
371 # missing if it's absent from all fields (as is the case when
372 # there's no data), and we restore the per-field missing cells in
373 # deserialized.
374 always_missing: set[CellIJ] = set(bounds.missing)
375 for field in aperture_correction_map.values():
376 always_missing &= field.bounds.missing
377 bounds = CellGridBounds(grid=bounds.grid, missing=frozenset(always_missing), bbox=bounds.bbox)
378 if any(field.unit is not None for field in aperture_correction_map.values()): 378 ↛ 379line 378 didn't jump to line 379 because the condition on line 378 was never true
379 raise ValueError("Aperture corrections should be dimensionless.")
380 table = astropy.table.Table(
381 rows=[cell_index.as_tuple() for cell_index in bounds.cell_indices()], names=["cell_i", "cell_j"]
382 )
383 good_cell_mask = np.ones(bounds.subgrid_size.as_tuple(), dtype=bool)
384 for cell_index in bounds.missing:
385 array_index = cell_index - bounds.subgrid_start
386 good_cell_mask[array_index.i, array_index.j] = False
387 for name, field in aperture_correction_map.items():
388 table.add_column(field._array[good_cell_mask], name=name, copy=False)
389 return CellApertureCorrectionMapSerializationModel(
390 table=archive.add_table(table, name="table"), bounds=bounds
391 )
393 def deserialize(self, archive: InputArchive[Any], **kwargs: Any) -> dict[str, CellField]:
394 if kwargs: 394 ↛ 395line 394 didn't jump to line 395 because the condition on line 394 was never true
395 raise InvalidParameterError(
396 f"Unrecognized parameters for cell aperture correction map: {set(kwargs.keys())}."
397 )
398 good_cell_mask = np.zeros(self.bounds.subgrid_size.as_tuple(), dtype=bool)
399 table = archive.get_table(self.table)
400 for tbl_ij, cell_index in zip(
401 table["cell_i", "cell_j"].iterrows(), self.bounds.cell_indices(), strict=True
402 ):
403 if cell_index.as_tuple() != tbl_ij: 403 ↛ 404line 403 didn't jump to line 404 because the condition on line 403 was never true
404 raise ArchiveReadError(
405 "Inconsistency between serialized aperture correction bounds and table."
406 )
407 array_index = cell_index - self.bounds.subgrid_start
408 good_cell_mask[array_index.i, array_index.j] = True
409 result: dict[str, CellField] = {}
410 for name, column in table.columns.items():
411 if name in ("cell_i", "cell_j"):
412 continue
413 extra_missing_mask = np.isnan(column)
414 # If there are any NaN entries in this column, make a custom
415 # CellGridBounds with additional missing cells for this field.
416 bounds = self.bounds
417 if extra_missing_mask.any():
418 extra_missing = frozenset(
419 {CellIJ(i=i, j=j) for i, j in table["cell_i", "cell_j"][extra_missing_mask].iterrows()}
420 )
421 bounds = CellGridBounds(
422 grid=self.bounds.grid, missing=self.bounds.missing | extra_missing, bbox=self.bounds.bbox
423 )
424 array = np.full(self.bounds.subgrid_size.as_tuple(), np.nan, dtype=np.float64)
425 array[good_cell_mask] = column
426 result[name] = CellField(bounds, array)
427 return result