Coverage for python/lsst/images/cells/_provenance.py: 68%
194 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-30 08:38 +0000
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__ = ("CoaddProvenance", "CoaddProvenanceSerializationModel")
16from collections.abc import Iterable
17from typing import TYPE_CHECKING, Any, ClassVar
19import astropy.table
20import astropy.units as u
21import numpy as np
22import pydantic
24from .._cell_grid import CellGridBounds, CellIJ
25from .._polygon import Polygon
26from ..describe import DescribableMixin, DescribeOptions, FieldRole, Report, ReportField
27from ..serialization import ArchiveTree, InputArchive, InvalidParameterError, OutputArchive, TableModel
29if TYPE_CHECKING:
30 try:
31 from lsst.afw.geom import Polygon as LegacyPolygon
32 from lsst.cell_coadds import CoaddInputs as LegacyCellCoaddInputs
33 from lsst.cell_coadds import MultipleCellCoadd as LegacyMultipleCellCoadd
34 from lsst.cell_coadds import ObservationIdentifiers as LegacyObservationIdentifiers
35 from lsst.skymap import Index2D as LegacyIndex2D
36 except ImportError:
37 type LegacyIndex2D = Any # type: ignore[no-redef]
38 type LegacyCellCoaddInputs = Any # type: ignore[no-redef]
39 type LegacyPolygon = Any # type: ignore[no-redef]
40 type LegacyMultipleCellCoadd = Any # type: ignore[no-redef]
41 type LegacyObservationIdentifiers = Any # type: ignore[no-redef]
44class CoaddProvenance(DescribableMixin):
45 """A pair of tables that record the inputs to a cell-based coadd.
47 Parameters
48 ----------
49 inputs
50 A table of {visit, detector} combinations that contribute to any cell
51 in the coadd.
52 contributions
53 A table of {visit, detector, cell} combinations that describe how an
54 observation contributed to a cell.
56 Notes
57 -----
58 This object can represent the provenance of a whole patch, a single cell,
59 or anything in between. In the single-cell case, the ``inputs`` and
60 ``contributions`` tables have the same number of rows (but may not be
61 ordered the same way!).
62 """
64 def __init__(self, inputs: astropy.table.Table, contributions: astropy.table.Table) -> None:
65 self._inputs = inputs
66 self._contributions = contributions
68 _INPUT_TABLE_COLUMNS: ClassVar[list[tuple[str, type, str]]] = [
69 ("instrument", np.object_, "Name of the instrument."),
70 ("visit", np.uint64, "ID of the visit."),
71 ("detector", np.uint16, "ID of the detector."),
72 ("physical_filter", np.object_, "Full name of the bandpass filter."),
73 ("day_obs", np.uint32, "Observation night as a YYYYMMDD integer."),
74 (
75 "polygon",
76 np.object_,
77 (
78 "Polygon that approximates the overlap of the observation and the coadd patch, "
79 "in coadd coordinates."
80 ),
81 ),
82 ]
84 _CONTRIBUTION_TABLE_COLUMNS: ClassVar[list[tuple[str, type, str, u.UnitBase | None]]] = [
85 ("cell_i", np.uint16, "Y-axis index of the cell within the patch.", None),
86 ("cell_j", np.uint16, "X-axis index of the cell within the patch.", None),
87 ("instrument", np.object_, "Name of the instrument.", None),
88 ("visit", np.uint64, "ID of the visit.", None),
89 ("detector", np.uint16, "ID of the detector.", None),
90 ("overlaps_center", np.bool_, "Whether a this observation overlaps the center of the cell.", None),
91 ("overlap_fraction", np.float64, "Fraction of the cell that is covered by the overlap region.", None),
92 ("unmasked_fraction", np.float64, "Fraction of the cell propagated to the coadd.", None),
93 ("weight", np.float64, "Weight to be used for this input in this cell.", None),
94 ("psf_shape_xx", np.float64, "Second order moments of the PSF.", u.pix**2),
95 ("psf_shape_yy", np.float64, "Second order moments of the PSF.", u.pix**2),
96 ("psf_shape_xy", np.float64, "Second order moments of the PSF.", u.pix**2),
97 (
98 "psf_shape_flag",
99 np.bool_,
100 "Flag indicating whether the PSF shape measurement was successful.",
101 None,
102 ),
103 ]
105 @classmethod
106 def make_empty_input_table(cls, n_rows: int) -> astropy.table.Table:
107 """Make an empty `inputs` table with a set number of rows.
109 Parameters
110 ----------
111 n_rows
112 Number of rows in the new table.
113 """
114 return astropy.table.Table(
115 [
116 astropy.table.Column(name=name, length=n_rows, dtype=dtype, description=description)
117 for name, dtype, description in cls._INPUT_TABLE_COLUMNS
118 ]
119 )
121 @classmethod
122 def make_empty_contribution_table(cls, n_rows: int) -> astropy.table.Table:
123 """Make an empty `contributions` table with a set number of rows.
125 Parameters
126 ----------
127 n_rows
128 Number of rows in the new table.
129 """
130 return astropy.table.Table(
131 [
132 astropy.table.Column(
133 name=name, length=n_rows, dtype=dtype, description=description, unit=unit
134 )
135 for name, dtype, description, unit in cls._CONTRIBUTION_TABLE_COLUMNS
136 ]
137 )
139 @property
140 def inputs(self) -> astropy.table.Table:
141 """A table of {visit, detector} combinations that contribute to any
142 cell in the coadd.
143 """
144 return self._inputs
146 @property
147 def contributions(self) -> astropy.table.Table:
148 """A table of {visit, detector, cell} combinations that describe how an
149 observation contributed to a cell.
150 """
151 return self._contributions
153 def __getitem__(self, cell: CellIJ) -> CoaddProvenance:
154 return self.subset([cell])
156 def subset(self, cells: Iterable[CellIJ]) -> CoaddProvenance:
157 """Return a new provenance object with just the given cells.
159 Parameters
160 ----------
161 cells
162 Cells to keep in the returned provenance.
163 """
164 cells_to_keep = astropy.table.Table(
165 rows=[(index.i, index.j) for index in cells],
166 names=["cell_i", "cell_j"],
167 dtype=[np.uint16, np.uint16],
168 )
169 contributions = astropy.table.join(self._contributions, cells_to_keep)
170 assert contributions.columns.keys() == {name for name, _, _, _ in self._CONTRIBUTION_TABLE_COLUMNS}
171 inputs = astropy.table.join(contributions["instrument", "visit", "detector"], self._inputs)
172 assert inputs.columns.keys() == {name for name, _, _ in self._INPUT_TABLE_COLUMNS}
173 return CoaddProvenance(inputs=inputs, contributions=contributions)
175 def _describe(
176 self,
177 options: DescribeOptions = DescribeOptions(),
178 /,
179 *,
180 bounds: CellGridBounds | None = None,
181 ) -> Report:
182 """Return a `Report` describing this provenance.
184 Parameters
185 ----------
186 options : `DescribeOptions`, optional
187 Rendering options. `DescribeOptions.brief` reports only the
188 number of input images, which no column scan is needed to count.
189 bounds : `CellGridBounds`, optional
190 Cells the image this provenance belongs to has data for. When
191 given, the number of cells with contributions is reported as a
192 fraction of them.
194 Notes
195 -----
196 The report summarizes the tables rather than rendering them. Rich
197 renders an embedded `astropy.table.Table` as plain text and never
198 consults the table's own ``_repr_html_``, so `inputs` and
199 `contributions` describe themselves better than this report could,
200 and a patch has tens of thousands of contribution rows.
201 """
202 n_inputs = len(self._inputs)
203 summary = (
204 f"CoaddProvenance({n_inputs} input image{'s' if n_inputs != 1 else ''})"
205 if n_inputs
206 else "CoaddProvenance(no input images)"
207 )
208 if options.brief:
209 return Report(type_name="CoaddProvenance", summary=summary)
210 fields: list[ReportField] = []
211 if n_inputs:
212 for column in ("instrument", "physical_filter"):
213 fields.append(
214 ReportField(
215 label=column,
216 value=", ".join(sorted({str(value) for value in self._inputs[column]})),
217 role=FieldRole.DERIVED,
218 )
219 )
220 n_visits = len(np.unique(self._inputs["visit"]))
221 input_images = f"{n_inputs} from {n_visits} visit{'s' if n_visits != 1 else ''}"
222 else:
223 input_images = "none"
224 fields.append(ReportField(label="input images", value=input_images, role=FieldRole.DERIVED))
225 if n_inputs:
226 first_night = int(self._inputs["day_obs"].min())
227 last_night = int(self._inputs["day_obs"].max())
228 fields.append(
229 ReportField(
230 label="day_obs",
231 value=str(first_night) if first_night == last_night else f"{first_night} - {last_night}",
232 role=FieldRole.DERIVED,
233 )
234 )
235 _, counts = np.unique(
236 np.column_stack([self._contributions["cell_i"], self._contributions["cell_j"]]),
237 axis=0,
238 return_counts=True,
239 )
240 if bounds is not None:
241 n_cells_with_data = bounds.subgrid_size.i * bounds.subgrid_size.j - len(bounds.missing)
242 cells = f"{len(counts)} of {n_cells_with_data} with contributions"
243 else:
244 cells = f"{len(counts)} with contributions" if len(counts) else "none"
245 fields.append(ReportField(label="cells", value=cells, role=FieldRole.DERIVED))
246 if len(counts):
247 low = int(counts.min())
248 high = int(counts.max())
249 per_cell = (
250 f"{low} input image{'s' if low != 1 else ''}"
251 if low == high
252 else f"{low} - {high} input images (median {np.median(counts):g})"
253 )
254 fields.append(ReportField(label="per cell", value=per_cell, role=FieldRole.DERIVED))
255 return Report(type_name="CoaddProvenance", summary=summary, fields=fields)
257 def serialize(self, archive: OutputArchive[Any]) -> CoaddProvenanceSerializationModel:
258 """Serialize the provenance to an output archive.
260 Parameters
261 ----------
262 archive
263 Archive to write to.
264 """
265 inputs = self._inputs.copy(copy_data=False)
266 contributions = self._contributions.copy(copy_data=False)
267 instrument = CoaddProvenanceSerializationModel._fix_str_for_serialization(
268 "instrument", inputs, contributions
269 )
270 physical_filter = CoaddProvenanceSerializationModel._fix_str_for_serialization(
271 "physical_filter", inputs
272 )
273 CoaddProvenanceSerializationModel._fix_polygon_for_serialization(inputs)
274 inputs_model = archive.add_table(inputs, name="inputs")
275 contributions_model = archive.add_table(contributions, name="contributions")
276 return CoaddProvenanceSerializationModel(
277 instrument=instrument,
278 physical_filter=physical_filter,
279 inputs=inputs_model,
280 contributions=contributions_model,
281 )
283 @staticmethod
284 def from_legacy(legacy_cell_coadd: LegacyMultipleCellCoadd) -> CoaddProvenance:
285 """Extract provenance from a legacy
286 `lsst.cell_coadds.MultipleCellCoadd` object.
288 Parameters
289 ----------
290 legacy_cell_coadd
291 Legacy cell coadd to extract provenance from.
292 """
293 inputs = CoaddProvenance.make_empty_input_table(len(legacy_cell_coadd.common.visit_polygons))
294 for n, (legacy_identifiers, legacy_polygon) in enumerate(
295 legacy_cell_coadd.common.visit_polygons.items()
296 ):
297 inputs["instrument"][n] = legacy_identifiers.instrument
298 inputs["visit"][n] = legacy_identifiers.visit
299 inputs["detector"][n] = legacy_identifiers.detector
300 inputs["physical_filter"][n] = legacy_identifiers.physical_filter
301 inputs["day_obs"][n] = legacy_identifiers.day_obs
302 inputs["polygon"][n] = Polygon.from_legacy(legacy_polygon)
303 n_contributions = 0
304 for legacy_cell in legacy_cell_coadd.cells.values():
305 n_contributions += len(legacy_cell.inputs)
306 contributions = CoaddProvenance.make_empty_contribution_table(n_contributions)
307 n = 0
308 for legacy_cell in legacy_cell_coadd.cells.values():
309 for legacy_identifiers, legacy_inputs in legacy_cell.inputs.items():
310 contributions["cell_i"][n] = legacy_cell.identifiers.cell.y
311 contributions["cell_j"][n] = legacy_cell.identifiers.cell.x
312 contributions["instrument"][n] = legacy_identifiers.instrument
313 contributions["visit"][n] = legacy_identifiers.visit
314 contributions["detector"][n] = legacy_identifiers.detector
315 contributions["overlaps_center"][n] = legacy_inputs.overlaps_center
316 contributions["overlap_fraction"][n] = legacy_inputs.overlap_fraction
317 contributions["unmasked_fraction"][n] = legacy_inputs.unmasked_overlap_fraction
318 contributions["weight"][n] = legacy_inputs.weight
319 contributions["psf_shape_xx"][n] = legacy_inputs.psf_shape.getIxx()
320 contributions["psf_shape_yy"][n] = legacy_inputs.psf_shape.getIyy()
321 contributions["psf_shape_xy"][n] = legacy_inputs.psf_shape.getIxy()
322 contributions["psf_shape_flag"][n] = legacy_inputs.psf_shape_flag
323 n += 1
324 return CoaddProvenance(inputs=inputs, contributions=contributions)
326 def to_legacy_polygon_map(self) -> dict[LegacyObservationIdentifiers, LegacyPolygon]:
327 """Construct a legacy mapping from
328 `lsst.cell_coadds.ObservationIdentifiers` to `lsst.afw.geom.Polygon`
329 from the `inputs` table.
330 """
331 from lsst.cell_coadds import ObservationIdentifiers as LegacyObservationIdentifiers
333 return {
334 LegacyObservationIdentifiers(
335 instrument=str(row["instrument"]),
336 physical_filter=str(row["physical_filter"]),
337 visit=int(row["visit"]),
338 day_obs=int(row["day_obs"]),
339 detector=int(row["detector"]),
340 ): row["polygon"].to_legacy()
341 for row in self.inputs
342 }
344 def to_legacy_cell_coadd_inputs(
345 self, observations: Iterable[LegacyObservationIdentifiers] | None
346 ) -> dict[LegacyIndex2D, dict[LegacyObservationIdentifiers, LegacyCellCoaddInputs]]:
347 """Construct a mapping from legacy cell index to the list of legacy
348 input structs for that cell.
350 Parameters
351 ----------
352 observations
353 Observations to include, or `None` to include all observations
354 in the `inputs` table.
355 """
356 from lsst.afw.geom.ellipses import Quadrupole
357 from lsst.cell_coadds import CoaddInputs as LegacyCoaddInputs
358 from lsst.skymap import Index2D as LegacyIndex2D
360 if observations is None:
361 observations = self.to_legacy_polygon_map().keys()
362 observations_by_key: dict[tuple[str, int, int], LegacyObservationIdentifiers] = {
363 (obs.instrument, obs.visit, obs.detector): obs for obs in observations
364 }
365 result: dict[LegacyIndex2D, dict[LegacyObservationIdentifiers, LegacyCoaddInputs]] = {}
366 for row in self.contributions:
367 obs_key = (str(row["instrument"]), int(row["visit"]), int(row["detector"]))
368 obs = observations_by_key[obs_key]
369 cell_inputs = result.setdefault(LegacyIndex2D(x=int(row["cell_j"]), y=int(row["cell_i"])), {})
370 cell_inputs[obs] = LegacyCoaddInputs(
371 overlaps_center=bool(row["overlaps_center"]),
372 overlap_fraction=float(row["overlap_fraction"]),
373 unmasked_overlap_fraction=float(row["unmasked_fraction"]),
374 weight=float(row["weight"]),
375 psf_shape=Quadrupole(row["psf_shape_xx"], row["psf_shape_yy"], row["psf_shape_xy"]),
376 psf_shape_flag=bool(row["psf_shape_flag"]),
377 )
378 return result
381class CoaddProvenanceSerializationModel(ArchiveTree):
382 """A Pydantic model used to represent a serialized `CoaddProvenance`.
384 Notes
385 -----
386 We can't rewrite the Astropy tables directly into the archive (e.g. as
387 FITS binary tables for a FITS archive), because:
389 - `str` columns are a huge pain in both Numpy and FITS;
390 - the polygon columns need to be rewritten as array-valued columns.
392 To deal with the string columns (``instrument`` and ``physical_filter``)
393 we do dictionary compression: we map each distinct value of those columns
394 to an integer, and then we save that mapping to the model while saving
395 an integer version of that column in the table. But if there is actually
396 only one value in that column (the most common case by far) we just drop
397 the column and store that value directly in the model.
398 """
400 SCHEMA_NAME: ClassVar[str] = "coadd_provenance"
401 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
402 MIN_READ_VERSION: ClassVar[int] = 1
403 PUBLIC_TYPE: ClassVar[type] = CoaddProvenance
405 instrument: str | dict[str, int] = pydantic.Field(
406 description=(
407 "Instrument name for all inputs to this coadd, or a mapping from "
408 "instrument name to the integer used in its place in the tables."
409 )
410 )
411 physical_filter: str | dict[str, int] = pydantic.Field(
412 description="Physical filter name for all inputs to this coadd."
413 )
414 inputs: TableModel = pydantic.Field(description="Table of all inputs to the coadd.")
415 contributions: TableModel = pydantic.Field(description="Table of per-cell contributions to the coadd.")
417 def deserialize(self, archive: InputArchive[Any], **kwargs: Any) -> CoaddProvenance:
418 """Deserialize a provenance from an input archive.
420 Parameters
421 ----------
422 archive
423 Archive to read from.
424 **kwargs
425 Unsupported keyword arguments are accepted only to provide
426 better error messages (raising
427 `.serialization.InvalidParameterError`).
429 Notes
430 -----
431 While `CoaddProvenance.subset` can be used to filter provenance
432 information down to just certain cells, there is no advantage to be
433 had from doing this during deserialization (the table data is not
434 ordered by cell, and hence there's read-slicing we can do).
435 """
436 if kwargs: 436 ↛ 437line 436 didn't jump to line 437 because the condition on line 436 was never true
437 raise InvalidParameterError(f"Unrecognized parameters for CoaddProvenance: {set(kwargs.keys())}.")
438 inputs = archive.get_table(self.inputs)
439 contributions = archive.get_table(self.contributions)
440 CoaddProvenanceSerializationModel._fix_str_for_deserialization(
441 "instrument", self.instrument, inputs, contributions
442 )
443 CoaddProvenanceSerializationModel._fix_str_for_deserialization(
444 "physical_filter", self.physical_filter, inputs
445 )
446 CoaddProvenanceSerializationModel._fix_polygon_for_deserialization(inputs)
447 for name, _, description in CoaddProvenance._INPUT_TABLE_COLUMNS:
448 inputs.columns[name].description = description
449 for name, _, description, unit in CoaddProvenance._CONTRIBUTION_TABLE_COLUMNS:
450 contributions.columns[name].description = description
451 contributions.columns[name].unit = unit
452 return CoaddProvenance(inputs=inputs, contributions=contributions)
454 @staticmethod
455 def _fix_str_for_serialization(column: str, *tables: astropy.table.Table) -> str | dict[str, int]:
456 """Rewrite a string column as an integer column or drop it.
458 Parameters
459 ----------
460 column
461 Name of the column to rewrite.
462 *tables
463 One or more astropy tables to rewrite. The first table is assumed
464 to have all values for this column that might appear in any other
465 tables.
467 Returns
468 -------
469 `str` | `dict` [`str`, `int`]
470 If there is only one unique value for this column in the first
471 table, that value (and the column will have been dropped from
472 all givne tables). If the tables are empty, the column is
473 dropped and an empty `dict` is returned. In all other cases the
474 given column is replaced with an integer column in all given
475 tables and the mapping from strings to integers is returned.
476 """
477 result: str | dict[str, int] = {name: n for n, name in enumerate(sorted(set(tables[0][column])))}
478 match len(result):
479 case 0: 479 ↛ 480line 479 didn't jump to line 480 because the pattern on line 479 never matched
480 pass
481 case 1: 481 ↛ 483line 481 didn't jump to line 483 because the pattern on line 481 always matched
482 (result,) = result.keys() # type: ignore[union-attr]
483 case _:
484 for table in tables:
485 table.columns[column] = astropy.table.Column(
486 data=[result[k] for k in table.columns[column]],
487 name=column,
488 dtype=np.uint8,
489 description=f"Integer mapped to {column} name.",
490 )
491 return result
492 # If we didn't remap to an integer (case 0 and 1 above), delete the
493 # column.
494 for table in tables:
495 del table.columns[column]
496 return result
498 @staticmethod
499 def _fix_str_for_deserialization(
500 column: str, value: str | dict[str, int], *tables: astropy.table.Table
501 ) -> None:
502 """Rewrite an integer column back to a string one.
504 Parameters
505 ----------
506 column
507 Name of the column to rewrite.
508 value
509 Value or mapping of values returned by
510 `_fix_str_for_serialization`.
511 tables
512 Tables to rewrite this column in.
513 """
514 match value:
515 case str(): 515 ↛ 518line 515 didn't jump to line 518 because the pattern on line 515 always matched
516 for table in tables:
517 table.columns[column] = astropy.table.Column([value] * len(table), dtype=object)
518 case dict():
519 mapping = {v: k for k, v in value.items()}
520 for table in tables:
521 table.columns[column] = astropy.table.Column(
522 [mapping[k] for k in table[column]], dtype=object
523 )
525 @staticmethod
526 def _fix_polygon_for_serialization(inputs: astropy.table.Table) -> None:
527 """Rewrite a polygon `object` column as a pair of array-valued columns
528 and an array-size column.
530 Parameters
531 ----------
532 inputs
533 A copy of the in-memory coadd inputs table to modify in-place into
534 its serialization form.
535 """
536 max_n_vertices = max(p.n_vertices for p in inputs["polygon"])
537 inputs["n_vertices"] = astropy.table.Column(
538 [p.n_vertices for p in inputs["polygon"]],
539 name="n_vertices",
540 dtype=np.uint8,
541 description="Number of polygon vertices.",
542 )
543 inputs["x_vertices"] = astropy.table.Column(
544 name="x_vertices",
545 dtype=np.float64,
546 length=len(inputs),
547 shape=(max_n_vertices,),
548 description="X coordinates of polygon vertices, in tract coordinates.",
549 )
550 inputs["x_vertices"][:, :] = np.nan
551 inputs["y_vertices"] = astropy.table.Column(
552 name="y_vertices",
553 dtype=np.float64,
554 length=len(inputs),
555 shape=(max_n_vertices,),
556 description="Y coordinates of polygon vertices, in tract coordinates.",
557 )
558 inputs["y_vertices"][:, :] = np.nan
559 for i, polygon in enumerate(inputs["polygon"]):
560 inputs["n_vertices"][i] = polygon.n_vertices
561 inputs["x_vertices"][i][: polygon.n_vertices] = polygon.x_vertices
562 inputs["y_vertices"][i][: polygon.n_vertices] = polygon.y_vertices
563 del inputs["polygon"]
565 @staticmethod
566 def _fix_polygon_for_deserialization(inputs: astropy.table.Table) -> None:
567 """Rewrite a a pair of array-valued columns and an array-size column
568 into a polygon `object` column.
570 Parameters
571 ----------
572 inputs
573 The serialized version of the coadd inputs table, to be modified
574 in-place into its in-memory form.
575 """
576 polygons = [
577 Polygon(x_vertices=x_vertices[:n_vertices], y_vertices=y_vertices[:n_vertices])
578 for n_vertices, x_vertices, y_vertices in zip(
579 inputs["n_vertices"], inputs["x_vertices"], inputs["y_vertices"]
580 )
581 ]
582 del inputs["n_vertices"]
583 del inputs["x_vertices"]
584 del inputs["y_vertices"]
585 inputs["polygon"] = astropy.table.Column(polygons, name="polygon", dtype=np.object_)